******************************************************************************************************
 /* ******************************************************
  Project:      Impact evaluation of the World Bank's Ethiopian Sustainable Land Management Project (SLMP) I (P090789) and II (P133133/P133410)
  
  Description:  This do file generates the ATT estimates 
  
  Created on: 01/24/22
  By: Susana Constenla 
 
 ************************************************************
 */
clear all
		mata: mata clear
		cap log close
		set more off
		set trace off
		set linesize 140
		set memory 200m
		set matsize 1000
		set maxvar 120000
		*set scheme plotplainifpri
		set seed 123456789
		timer clear
		label drop _all 
		estimates clear
graph set window fontface "Times New Roman"

********************************************************************************
* DEFINITION OF THE WORKING DIRECTORY. INPUT AND OUTPUT FILES*

		if "`c(os)'"=="Windows" {								
		local drive: env HOMEDRIVE
		local path : env HOMEPATH
			global 	BASE_DIR "`drive'`path'\PATH_TO_MAIN_FOLDER" // change for path to main working directory
			}
		if "`c(os)'"=="MacOSX" {
		local path : env HOMEPATH
			global BASE_DIR "~//`path'\PATH_TO_MAIN_FOLDER" // change for path to main working directory
		}

*----------------------------
cd "$BASE_DIR"
 
 *IN*
		global INPUT	"$BASE_DIR\INPUT_DATA_PATH" //change for input data path

 *OUT*	
        global GRAPHS   "$BASE_DIR\OUTPUT_GRAPH_PATH"  //change for output graph path
	    global TABLES   "$BASE_DIR\OUTPUT_TABLES_PATH" //change for output tables path
		
****************************************************************************************************************

	use "$INPUT/SLMP_EVI_cleaned_pixel_level.dta", replace 

* ############## Creating new vars to do the analysis by season ##################### * 

	tostring year, gen(year_s) 
	tostring month, gen(month_s)
	
gen treated=(donor=="WBI")
lab define treated 1 "Treated[SLMP1]" 0 "Control[SLMP2]"
la val treated treated 

 
* Creating season variables to average satelite indicators over - following FAO's seasonality classification for Ethiopia
 gen season=.
	replace season=1 if (month>=10 | month==1)
	replace season=2 if  month>=2 & month<6
	replace season=3 if  month>=6 & month<10
lab define season 1 "Dry Season" 2 "Pre-rainy season" 3 "Rainy season"

la val season season 
  tostring season, gen(season_s)
  gen year_season_s=year_s+season_s
  destring year_season_s, gen(year_season)
 
 rename evi_cleaned evi_pixel 

 *-----------------------------------------*
 * Averaging variables at the season level * 
 *-----------------------------------------*
 local vars  total_precip_chirpsv20 av_temp_era5 
 foreach var of varlist `vars' {
  bys mws_pixel_id_1km year_season: egen `var'_s=mean(`var') // averaging at the year_season level
  }
  
  * Unbalanced panel due to cleaning: Just take into account pixels for which there is a value for each month of the season *
 preserve 
	  sum evi_pixel
	   gen nm=cond(evi_pixel!=.,0,1)
	   bys mws_pixel_id_1km year_season: egen missing=total(nm)
	   drop if missing!=0
	   
	  foreach var of varlist evi_pixel {
	  bys mws_pixel_id_1km year_season: egen evi_pixel_s=mean(evi_pixel) // averaging at the year_season level
	  }
	   bys mws_pixel_id_1km year_season: egen evi_max_s=max(evi_pixel) // max of the the year_season level 

	  keep mws_pixel_id_1km year_season evi_pixel_s  evi_max_s 
	  bys mws_pixel_id_1km year_season: keep if _n==1
	  tempfile average
	  save `average', replace 
 restore 
   
  merge m:1 mws_pixel_id_1km year_season using `average', nogen
  
 * Checking out missing values on evi after cleaning, by group *
  bys mws_pixel_id_1km: egen non_missing=total(evi_pixel!=.)
  replace non_missing=(non_missing/248)*100
  bys treated: sum non_missing if season==3, de
  
  ******************************
   drop if mws_pixel_id_1km==. //76 obs for which there is not EVI data (not spatial join)

   
    gen crop_range_mask=(crop_mask_v2+rangeland_mask_v2)
    gen crop_range_frac=(crop_frac+range_frac)
	
	sum crop* range*
*---------------------------------------------------------------------------------
   
   gen drought_SPI_1=(SPI<=-1) if SPI!=.
   gen drought_SPI_15=(SPI<=-1.5) if SPI!=.
   gen drought_SPI_2=(SPI<=-2) if SPI!=.

    foreach var of varlist drought* {
   bys mws_pixel_id_1km year_season: egen max_`var'=max(`var')
   replace `var'= max_`var'
   drop max_`var'
  }
  
********************************************************
* GRAPH of drought over time *
  	   preserve
		keep if elevation>1000 & (crop_mask_v2>50|rangeland_mask_v2>50|crop_range_mask>50)
		bys mws_pixel_id_1km year_season: keep if _n==1
		
		foreach var of varlist drought_SPI_1 drought_SPI_15 drought_SPI_2{
			bys year_season: egen prop_`var'_1 = mean(`var') if treated==1
			bys year_season: egen prop_`var'_0 = mean(`var') if treated==0
		
		
  twoway (line  prop_`var'_1  year if treated==1 & season==3 & year>=2002 & year<=2015, lcolor(orange red) lwidth(medthick) ///
         ylabel(0(0.2)1, format(%03.2f) labsize(medium) glpattern(dot) glcolor(gray%90) glwidth(thick) ) ///
		 legend(lab(1 "SLMP1") lab(2 "SLMP2") col(2) position(6) size(medium)) ///
		 xlabel(#16,  glpattern(dot) glcolor(gray%90) glwidth(thick) labsize(medium)) xtitle("", height(10) ) ///
		 ytitle("Prop. of pixels that suffered a drought event", height(10) size(medlarge)) ) ///
		 (line  prop_`var'_0 year if treated==0 & season==3 & year>=2002 & year<=2015, lcolor(black%50) lpattern(dash) lwidth(medthick))
  graph export "$GRAPHS/lpolyreg_`var'_1km_season3.tif", replace
		}
     restore 
	 
 *============================================================================================
 *##### UNBALANCED PANEL GRAPH  ######*
 *------------------------------------
 preserve
  keep if season==3 & evi_pixel_s!=. //just pixels that form part of the estimation model
  bys mws_pixel_id_1km year_season: keep if _n==1
  keep if year<=2015 & year>=2002
 
*--------------------------

 unique mws_pixel_id_1km if elevation>1000 & (crop_mask_v2>50|rangeland_mask_v2>50|crop_range_mask>50), by(treated) gen(total_pixels) 
   tab total_pixels

	gen prop_pixels_0=.
	gen prop_pixels_1=.
	
	 levelsof year_season, local(year_season)

 foreach l of local year_season {
  count  if treated==1 & year_season==`l' & elevation>1000 & (crop_mask_v2>50|rangeland_mask_v2>50|crop_range_mask>50)
   replace prop_pixels_1=`r(N)'/2052 if year_season==`l' & treated==1
   count  if treated==0 & year_season==`l' &  elevation>1000 & (crop_mask_v2>50|rangeland_mask_v2>50|crop_range_mask>50) 
   replace prop_pixels_0=`r(N)'/4392  if year_season==`l' & treated==0
   
    }
	replace prop_pixels_0=prop_pixels_0*100
 	replace prop_pixels_1=prop_pixels_1*100

	sort year 
	twoway (line prop_pixels_1 year, lwidth(medthick) lcolor(orange red)) (line prop_pixels_0 year,  lcolor(black) lwidth(medthick) ///
	       xlabel(#16,labsize(medium) glpattern(dot) glcolor(gray%90) glwidth(thick)) ///
		   legend(lab(1 "SLMP1") lab(2 "SLMP2") col(2) position(5) size(medium)) ///
		   xtitle("", height(5)) ytitle("% of total pixels with valid EVI data", size(medlarge)  height(8)) ///
		   ylabel(0(10)70, labsize(medium) glpattern(dot) glcolor(gray%90) glwidth(thick)) ///
		   note("Total number of distinct pixels included in the panel:    SLMP1:2,052    SLMP2: 4,392", size(small)))
	
	   graph export "$GRAPHS/SI_Figure6_Panel_balance.tif", replace

	*---------------------------------
restore 
	
*--------------------------
* NOTE: sum stats of time invariant vars (soil quality, elevation and crop mask) is in another dofile, with the 3-group balance comparison
*---------------------------------
 *=========================================================================
 * APPROACH 1: PIPELINE DESIGN, SLMP1 = TREATMENT ; SLMP2 = CONTROL
 *=========================================================================

 *---------------------------------------------
 * 1. Cheking for paralell trends assumption 
 *--------------------------------------------
		*1.1. Visual approach - local polynomial non parametric regression of EVI over time by treatment group
  
  *---------- 
  preserve
	  keep if season==3 & year<=2015 & year>=2002

	  bys year_season: gen time_cons=1 if _n==1 
	  replace time_cons=sum(time_cons)
	  replace time_cons=time_cons-1

	
	  twoway (lpolyci  evi_pixel_s time_cons if treated==1 & elevation>1000 & (crop_mask_v2>50|rangeland_mask_v2>50|crop_range_mask>50) , bwidth(0.5)  ///
	legend(order(2 "SLMP1" 4 "SLMP2") col(2)  position(6) height(5) size(medium))   xline(8.5 13) ///
	ytitle("EVI", height(8) size(medlarge)) lwidth(medthick) fcolor(orange red) fintensity(inten10) ///
	xlabel(0 "2002" 1 "2003" 2 "2004" 3 "2005" 4 "2006" 5 "2007" 6 "2008" 7 "2009" 8 "2010" 9 "2011" 10 "2012" 11 "2013" 12 "2014" 13 "2015",  glpattern(dot) glcolor(gray%90) glwidth(thick) labsize(medium)) ///
	ylabel(0.28(0.02)0.40, glpattern(dot) glcolor(gray%90) glwidth(thick) format(%03.2f)   labsize(medium)))   ///
	(lpolyci  evi_pixel_s time_cons if treated==0 & elevation>1000 & (crop_mask_v2>50|rangeland_mask_v2>50|crop_range_mask>50)  ,  bwidth(0.5)   ///
	  lcolor(black)  alcolor(black%30) ciplot(rline) alpattern(solid) fcolor(orange red) fintensity(inten10) lpattern(shortdash) xtitle("") lwidth(medthick))
	 graph export "$GRAPHS/Fig2_leftpanel_EVI_lpolyreg_season3.tif", replace
  
  restore 
  
 *---------------------------------------------------------
 
    gen post=(year>2010) // based on budget disbursement 
	gen post1=(year>2010 & year<=2013)
	gen post2=(year>=2014 & year<=2015)

	
 *=========================================================================
 *  SPECIFICATION for parallel trends check : EVENT STUDY
 *=========================================================================
  
   bys mws_pixel_id_1km year_season: keep if _n==1
   drop if mws_pixel_id_1km==. 
   
     	   bysort mws_pixel_id_1km (year_season) : gen total_precip_prev_s = (total_precip_chirpsv20_s[_n-1])

	  gen elev_preci=total_precip_chirpsv20_s*elevation
		 gen elev_temp=av_temp_era5_s*elevation
		  gen temp_sq = av_temp_era5_s^2
		  gen precip_sq = total_precip_chirpsv20_s^2
		  gen precip_wate= total_precip_chirpsv20_s*waterholdi
	  
	  
	  foreach var of varlist af* {
	  gen prec_`var'=`var'*total_precip_chirpsv20_s
	  gen temp_`var'=`var'*av_temp_era5_s	  
	  }

 global soil "prec_af250m_nut prec_af250m_n_1 prec_af250m_n_2  prec_af_pH515cm prec_af_ORCDRC1 prec_af_EC515cm  prec_af_Depthto prec_af_CEC515c precip_wate"

*#######################################################

tab year_season if year<=2015 & season==3, gen(period)
 foreach var of varlist period*{
  replace `var'=`var'*treated 
 }
 
 xtset mws_pixel_id_1km year_season
 
 la var period1 "period -11 (year 2000)"
 la var period2 "period -10  (year 2001)"
 la var period3 "period -9  (year 2002)"
 la var period4 "period -8  (year 2003)"
 la var period5 "period -7  (year 2004)"
 la var period6 "period -6  (year 2005)"
 la var period7 "period -5  (year 2006)"
 la var period8 "period -4  (year 2007)"
 la var period9 "period -3  (year 2008)"
 la var period10 "period -2  (year 2009)"
 la var period12 "period +1  (year 2011)"
 la var period13 "period +2  (year 2012)"
 la var period14 "period +3  (year 2013)"
 la var period15 "period +4  (year 2014)"
 la var period16 "period +5  (year 2015)"


	 constraint define 1   period3 + period4 + period5 ///
 + period6 + period7 + period8 + period9 + period10 + period11= 0
	

foreach var of varlist   evi_pixel_s  {

	cnsreg `var' ///
       period3-period16 ///
       i.mws_pixel_id_1km i.year_season ///
			 if  year<=2015 & year>=2002 & season==3 & elevation>1000 & (crop_mask_v2>50|rangeland_mask_v2>50|crop_range_mask>50), noconstant constraint(1) vce(cluster mws_pixel_id_1km) 
  outreg2 using "$TABLES/Event_Study_results_`var'_Millermethod.tex", keep(period*) replace label
  
  predict fitted if e(sample)
  corr `var' fitted if e(sample)
  drop fitted
 *--------------------------------------------------------------------------- 
*----------------------------------------------------------------------------  
  cnsreg `var' ///
       period3-period16 ///
	     total_precip_chirpsv20_s av_temp_era5_s ///
		 total_precip_prev_s elev_temp elev_preci temp_sq precip_sq  /// 
		        i.mws_pixel_id_1km i.year_season ///
			 if  year<=2015 & year>=2002 & season==3 & elevation>1000 & (crop_mask_v2>50|rangeland_mask_v2>50|crop_range_mask>50), noconstant constraint(1) vce(cluster mws_pixel_id_1km) 
  outreg2 using "$TABLES/Event_Study_results_`var'_Millermethod.tex", keep(period*) append label
  
    predict fitted if e(sample)
  corr `var' fitted if e(sample)
  drop fitted
 *----------------------------------------------------------------------------
  cnsreg `var' ///
       period3-period16 ///
	     total_precip_chirpsv20_s av_temp_era5_s ///
		 total_precip_prev_s elev_temp elev_preci temp_sq precip_sq  /// 
		 $soil ///
		 i.mws_pixel_id_1km i.year_season ///
			 if  year<=2015 & year>=2002 & season==3 & elevation>1000 & (crop_mask_v2>50|rangeland_mask_v2>50|crop_range_mask>50), noconstant constraint(1) vce(cluster mws_pixel_id_1km) 
  outreg2 using "$TABLES/Event_Study_results_`var'_Millermethod.tex", keep(period*) append label
  
    predict fitted if e(sample)
  corr `var' fitted if e(sample)
  drop fitted
 *----------------------------------------------------------------------------
  cnsreg `var' ///
       period3-period16 ///
	     total_precip_chirpsv20_s av_temp_era5_s ///
		 total_precip_prev_s elev_temp elev_preci temp_sq precip_sq  /// 
		 $soil ///
		 i.mws_pixel_id_1km i.year_season##i.aez_09 ///
			 if  year<=2015 & year>=2002 & season==3 & elevation>1000 & (crop_mask_v2>50|rangeland_mask_v2>50|crop_range_mask>50), noconstant constraint(1) vce(cluster mws_pixel_id_1km) 
			 
  outreg2 using "$TABLES/Event_Study_results_`var'_Millermethod.tex", keep(period*) append label
  
 *----------------------------------------------------------------------------

   foreach i in  3 4 5 6 7 8 9 10 11  12 13 14 15 16 {
    gen beta_`i'=_b[period`i'] if period`i'==1
	gen se_`i'=_se[period`i'] if period`i'==1
   }
   
   egen beta=rowmax(beta_*)
   egen se=rowmax(se_*)

   sort year_season

   serrbar beta se year if year<=2015 & year>=2002, scale (1.96)  xlabel(#15) ///
   ytitle("EVI", height(8) size(medlarge)) xtitle("") ///
   yline(0, lwidth(medthick)) mvopts(c(i)  m(D) msize(medium) mlw(medlarge) mlc(black%20) mc(gray%50)) ///
   xlabel(, glpattern(dot) glcolor(gray%90) glwidth(thick)  labsize(medium)) ylabel(, glpattern(dot) glcolor(gray%90) glwidth(thick)  format(%03.2f) labsize(medium)) ///
    lwidth(medthick) lcolor(orange_red) xline(2010.50, lcolor(orange_red) lwidth(medthick) lpattern(longdash)) 
  
       graph export "$GRAPHS/Fig3_panelA_evi_Event_Study_Millermethod.tif", replace
		
	drop beta* se se_* 

}


*================================================
* Robustness check - Alternative land cover map *

foreach var of varlist   evi_pixel_s  {

	cnsreg `var' ///
       period3-period16 ///
       i.mws_pixel_id_1km i.year_season ///
			 if  year<=2015 & year>=2002 & season==3 & elevation>1000 & (crop_frac>0.50|range_frac>0.50|crop_range_frac>0.50), noconstant constraint(1) vce(cluster mws_pixel_id_1km) 
  outreg2 using "$TABLES/Event_Study_results_`var'_altcover.tex", keep(period*) replace label
  
  predict fitted if e(sample)
  corr `var' fitted if e(sample)
  drop fitted
 *--------------------------------------------------------------------------- 
*----------------------------------------------------------------------------  
  cnsreg `var' ///
       period3-period16 ///
	     total_precip_chirpsv20_s av_temp_era5_s ///
		 total_precip_prev_s elev_temp elev_preci temp_sq precip_sq  /// 
		        i.mws_pixel_id_1km i.year_season ///
			 if  year<=2015 & year>=2002 & season==3 & elevation>1000 & (crop_frac>0.50|range_frac>0.50|crop_range_frac>0.50), noconstant constraint(1) vce(cluster mws_pixel_id_1km) 
  outreg2 using "$TABLES/Event_Study_results_`var'_altcover.tex", keep(period*) append label
  
    predict fitted if e(sample)
  corr `var' fitted if e(sample)
  drop fitted
 *----------------------------------------------------------------------------
  cnsreg `var' ///
       period3-period16 ///
	     total_precip_chirpsv20_s av_temp_era5_s ///
		 total_precip_prev_s elev_temp elev_preci temp_sq precip_sq  /// 
		 $soil ///
		 i.mws_pixel_id_1km i.year_season ///
			 if  year<=2015 & year>=2002 & season==3 & elevation>1000 & (crop_frac>0.50|range_frac>0.50|crop_range_frac>0.50), noconstant constraint(1) vce(cluster mws_pixel_id_1km) 
  outreg2 using "$TABLES/Event_Study_results_`var'_altcover.tex", keep(period*) append label
  
    predict fitted if e(sample)
  corr `var' fitted if e(sample)
  drop fitted
 *----------------------------------------------------------------------------
  cnsreg `var' ///
       period3-period16 ///
	     total_precip_chirpsv20_s av_temp_era5_s ///
		 total_precip_prev_s elev_temp elev_preci temp_sq precip_sq  /// 
		 $soil ///
		 i.mws_pixel_id_1km i.year_season##i.aez_09 ///
			 if  year<=2015 & year>=2002 & season==3 & elevation>1000 & (crop_frac>0.50|range_frac>0.50|crop_range_frac>0.50), noconstant constraint(1) vce(cluster mws_pixel_id_1km) 
			 
  outreg2 using "$TABLES/Event_Study_results_`var'_altcover.tex", keep(period*) append label
  
 *----------------------------------------------------------------------------

   foreach i in  3 4 5 6 7 8 9 10 11  12 13 14 15 16 {
    gen beta_`i'=_b[period`i'] if period`i'==1
	gen se_`i'=_se[period`i'] if period`i'==1
   }
   
   egen beta=rowmax(beta_*)
   egen se=rowmax(se_*)

   sort year_season

   serrbar beta se year if year<=2015 & year>=2002, scale (1.96)  xlabel(#15) ///
   ytitle("EVI", height(8) size(medlarge)) xtitle("") ///
   yline(0, lwidth(medthick)) mvopts(c(i)  m(D) msize(medium) mlw(medlarge) mlc(black%20) mc(gray%50)) ///
   xlabel(, glpattern(dot) glcolor(gray%90) glwidth(thick)  labsize(medium)) ylabel(, glpattern(dot) glcolor(gray%90) glwidth(thick)  format(%03.2f) labsize(medium)) ///
    lwidth(medthick) lcolor(orange_red) xline(2010.50, lcolor(orange_red) lwidth(medthick) lpattern(longdash)) 
  
       graph export "$GRAPHS/SI_Figure4_panelA_EVI_Event_Study_altcover.tif", replace
		
	drop beta* se se_* 

}

*=============================================================================================
  *=========================================================================
 *  DIFFERENCE IN DIFFERENCE
 *=========================================================================


 	  drop if year<2002 // to make equal to SIF series that are just available from this year on

 forvalues s=3/3{
 foreach var of varlist  evi_pixel_s  {
	  	  
	   xtset mws_pixel_id_1km year_season  
			
  *###--------------------------------
  
xtreg `var' ///
       treated##post i.year_season ///
			 if  year<=2015 & season==`s' & elevation>1000 &  (crop_mask_v2>50|rangeland_mask_v2>50|crop_range_mask>50), fe vce(cluster mws_pixel_id_1km)
			*------------------------------------------------
   outreg2 using "$TABLES/DID_results_`var'_season`s'.tex", ///
   keep(1.treated#1.post total_precip_chirpsv20_s av_temp_era5_s total_precip_prev_s elev_temp temp_sq precip_sq elev_preci) replace label
   
xtreg `var' ///
       treated##post i.year_season total_precip_chirpsv20_s av_temp_era5_s total_precip_prev_s ///
							 elev_temp elev_preci temp_sq precip_sq  ///
			 if  year<=2015 & season==`s' & elevation>1000 &  (crop_mask_v2>50|rangeland_mask_v2>50|crop_range_mask>50), fe vce(cluster mws_pixel_id_1km)
			*------------------------------------------------
   outreg2 using "$TABLES/DID_results_`var'_season`s'.tex", ///
   keep(1.treated#1.post total_precip_chirpsv20_s av_temp_era5_s total_precip_prev_s elev_temp temp_sq precip_sq elev_preci) append label

xtreg `var' ///
       treated##post i.year_season total_precip_chirpsv20_s av_temp_era5_s total_precip_prev_s ///
							 elev_temp elev_preci temp_sq precip_sq  ///
							 $soil ///
			 if  year<=2015 & season==`s' & elevation>1000 &  (crop_mask_v2>50|rangeland_mask_v2>50|crop_range_mask>50), fe vce(cluster mws_pixel_id_1km)
			*------------------------------------------------
   outreg2 using "$TABLES/DID_results_`var'_season`s'.tex", ///
   keep(1.treated#1.post total_precip_chirpsv20_s av_temp_era5_s total_precip_prev_s elev_temp temp_sq precip_sq elev_preci) append label
	
 xtreg `var' ///
       treated##post i.year_season##i.aez_09 total_precip_chirpsv20_s av_temp_era5_s total_precip_prev_s ///
							 elev_temp elev_preci temp_sq precip_sq  ///
							 $soil ///
			 if  year<=2015 & season==`s' & elevation>1000 &  (crop_mask_v2>50|rangeland_mask_v2>50|crop_range_mask>50), fe vce(cluster mws_pixel_id_1km)
			*------------------------------------------------
   outreg2 using "$TABLES/DID_results_`var'_season`s'.tex", ///
   keep(1.treated#1.post total_precip_chirpsv20_s av_temp_era5_s total_precip_prev_s elev_temp temp_sq precip_sq elev_preci) append label addnote(treatment group mean in post period == `m')
 
*===============================================================
* Coefficient storage for plotting it later *

 sum evi_pixel_s if treated==1 & post ==0 & e(sample)==1
 local m = `r(mean)'
 global scale1_EVI_`s' = (1/(`m'))*100 //% change formula
 
 mat evi_general_season`s' = _b[1.treated#1.post]
 
 local lb_post = _b[1.treated#1.post] - invttail(e(df_r),0.025)*_se[1.treated#1.post] 
 local ub_post = _b[1.treated#1.post] + invttail(e(df_r),0.025)*_se[1.treated#1.post]
 
 mat evi_general_ci_season`s' =  ((`lb_post')\(`ub_post')) 
 

*=====================================
* Dynamic program effects *

xtreg `var' ///
       treated##post1 treated##post2 i.year_season ///
			 if  year<=2015 & season==`s' & elevation>1000 &  (crop_mask_v2>50|rangeland_mask_v2>50|crop_range_mask>50), fe vce(cluster mws_pixel_id_1km)
			*------------------------------------------------
  outreg2 using "$TABLES/DID_results_`var'_2phases_season`s'.tex", ///
  keep(1.treated#1.post1 1.treated#1.post2 total_precip_chirpsv20_s av_temp_era5_s total_precip_prev_s elev_temp temp_sq precip_sq elev_preci) replace label

xtreg `var' ///
       treated##post1 treated##post2 i.year_season total_precip_chirpsv20_s av_temp_era5_s total_precip_prev_s ///
							 elev_temp elev_preci temp_sq precip_sq  ///
			 if  year<=2015 & season==`s' & elevation>1000 &  (crop_mask_v2>50|rangeland_mask_v2>50|crop_range_mask>50), fe vce(cluster mws_pixel_id_1km)
			*------------------------------------------------
  outreg2 using "$TABLES/DID_results_`var'_2phases_season`s'.tex", ///
  keep(1.treated#1.post1 1.treated#1.post2 total_precip_chirpsv20_s av_temp_era5_s total_precip_prev_s elev_temp temp_sq precip_sq elev_preci) append label

xtreg `var' ///
       treated##post1 treated##post2 i.year_season total_precip_chirpsv20_s av_temp_era5_s total_precip_prev_s ///
							 elev_temp elev_preci temp_sq precip_sq  ///
							 $soil ///
			 if  year<=2015 & season==`s' & elevation>1000 &  (crop_mask_v2>50|rangeland_mask_v2>50|crop_range_mask>50), fe vce(cluster mws_pixel_id_1km)
			*------------------------------------------------
  outreg2 using "$TABLES/DID_results_`var'_2phases_season`s'.tex", ///
  keep(1.treated#1.post1 1.treated#1.post2 total_precip_chirpsv20_s av_temp_era5_s total_precip_prev_s elev_temp temp_sq precip_sq elev_preci) append label

  
xtreg `var' ///
       treated##post1 treated##post2 i.year_season##i.aez_09 total_precip_chirpsv20_s av_temp_era5_s total_precip_prev_s ///
							 elev_temp elev_preci temp_sq precip_sq  ///
							 $soil ///
			 if  year<=2015 & season==`s' & elevation>1000 &  (crop_mask_v2>50|rangeland_mask_v2>50|crop_range_mask>50), fe vce(cluster mws_pixel_id_1km)
			*------------------------------------------------
	    outreg2 using "$TABLES/DID_results_`var'_2phases_season`s'.tex", ///
  keep(1.treated#1.post1 1.treated#1.post2 total_precip_chirpsv20_s av_temp_era5_s total_precip_prev_s elev_temp temp_sq precip_sq elev_preci) append label 
	
	estimates store DID_EVI_2

	*--------------------------------------------------------------------------

 global scale2_EVI_1_`s' = (1/(`m'))*100 //% change formula
 global scale2_EVI_2_`s'= (1/(`m'))*100 //% change formula
 
 mat evi_general_season`s' = nullmat(evi_general_season`s'), (_b[1.treated#1.post1]),  (_b[1.treated#1.post2])
 mat colnames evi_general_season`s' = ATE ATE_post1 ATE_post2
 
 forvalues i=1/2{
 local lb_post`i' = _b[1.treated#1.post`i'] - invttail(e(df_r),0.025)*_se[1.treated#1.post`i'] 
 local ub_post`i' = _b[1.treated#1.post`i'] + invttail(e(df_r),0.025)*_se[1.treated#1.post`i']
 }
 
 mat evi_general_ci_season`s' = nullmat(evi_general_ci_season`s'), ((`lb_post1')\(`ub_post1')) , ((`lb_post2')\(`ub_post2'))
  mat colnames evi_general_ci_season`s' = ATE ATE_post1 ATE_post2

  test _b[1.treated#1.post1]=_b[1.treated#1.post2]
  test _b[1.treated#1.post1]=0 ,notest 
  test _b[1.treated#1.post2]=0 , accum

*----------------------------------
 
*#################################################################################################################

 *=========================================================================
 *  DROUGHT ANALYSIS 
 *=========================================================================
 	   xtset mws_pixel_id_1km year_season 

	foreach x in   drought_SPI_15  drought_SPI_2   {
	xtreg `var' ///
       treated##`x'##post treated##`x' i.year_season ///
	if  year<=2015 & season==`s' & elevation>1000 &  (crop_mask_v2>50|rangeland_mask_v2>50|crop_range_mask>50), fe vce(cluster mws_pixel_id_1km)
	
	 outreg2 using "$TABLES/DID_results_`x'_`var'_season`s'.tex",keep(1.treated#1.`x'#1.post 1.treated#1.post 1.`x'#1.post 1.treated#1.`x' 1.`x' av_temp_era5_s  elev_temp temp_sq ) replace label
*-----------------------------------------------
 	xtreg `var' ///
       treated##`x'##post treated##`x' i.year_season   ///
	   	   	   	    av_temp_era5_s elev_temp temp_sq  ///
	if  year<=2015 & season==`s' & elevation>1000 &  (crop_mask_v2>50|rangeland_mask_v2>50|crop_range_mask>50), fe vce(cluster mws_pixel_id_1km)
	
	 outreg2 using "$TABLES/DID_results_`x'_`var'_season`s'.tex",keep(1.treated#1.`x'#1.post 1.treated#1.post 1.`x'#1.post 1.treated#1.`x' 1.`x' av_temp_era5_s  elev_temp temp_sq ) append label
*-----------------------------------------------
 	xtreg `var' ///
       treated##`x'##post treated##`x' i.year_season   ///
	   	   	   	    av_temp_era5_s elev_temp temp_sq  ///
					$soil ///
	if  year<=2015 & season==`s' & elevation>1000 &  (crop_mask_v2>50|rangeland_mask_v2>50|crop_range_mask>50), fe vce(cluster mws_pixel_id_1km)
	
	 outreg2 using "$TABLES/DID_results_`x'_`var'_season`s'.tex",keep(1.treated#1.`x'#1.post 1.treated#1.post 1.`x'#1.post 1.treated#1.`x' 1.`x' av_temp_era5_s  elev_temp temp_sq ) append label
*-----------------------------------------------
 	xtreg `var' ///
       treated##`x'##post treated##`x' i.year_season##i.aez_09   ///
	   	   	   	    av_temp_era5_s elev_temp temp_sq  ///
					$soil ///
	if  year<=2015 & season==`s' & elevation>1000 &  (crop_mask_v2>50|rangeland_mask_v2>50|crop_range_mask>50), fe vce(cluster mws_pixel_id_1km)
	  	 
	outreg2 using "$TABLES/DID_results_`x'_`var'_season`s'.tex",keep(1.treated#1.`x'#1.post 1.treated#1.post 1.`x'#1.post 1.treated#1.`x' 1.`x' av_temp_era5_s  elev_temp temp_sq ) append label
	
*============
	
	 sum evi_pixel_s if treated==1 & post ==0 & e(sample)==1
 local m1 = `r(mean)'
 
  global s_`x'_post_EVI = (1/(`m1'))*100 //% change formula
  global s_`x'_EVI = (1/(`m1'))*100 //% change formula
  global s_no_`x'_EVI = (1/(`m1'))*100 //% change formula

  
  lincom 1.treated#1.`x'#1.post+1.treated#1.post
   mat EVI_d_post=`r(estimate)', `r(se)', `r(p)', `r(lb)', `r(ub)'
   mat evi_`x' = `r(estimate)'
   mat evi_`x'_ci =  (`r(lb)' \ `r(ub)')

  mat evi_no_drought = (_b[1.treated#1.post])

 local lb_post = _b[1.treated#1.post] - invttail(e(df_r),0.025)*_se[1.treated#1.post] 
 local ub_post = _b[1.treated#1.post] + invttail(e(df_r),0.025)*_se[1.treated#1.post]
 
 mat evi_no_drought_ci = ((`lb_post')\(`ub_post')) 
 
  
*-----------------------------------------------	
*=======================================================

	 xtreg `var' ///
        treated##post1 treated##`x'##post1 treated##`x'##post2 treated##`x' i.year_season ///
	if  year<=2015 & season==`s' & elevation>1000 &  (crop_mask_v2>50|rangeland_mask_v2>50|crop_range_mask>50), fe vce(cluster mws_pixel_id_1km)

   outreg2 using "$TABLES/DID_results_`x'_`var'_2phases_season`s'.tex", ///
   keep(1.treated#1.post1 1.treated#1.post2 1.treated#1.`x'#1.post1 1.treated#1.`x'#1.post2 1.post1 1.post2  1.`x'#1.post1 1.`x'#1.post2 1.treated#1.`x' 1.`x' av_temp_era5_s  elev_temp temp_sq ) append label
*-----------------------------------------------   
   	 xtreg `var' ///
        treated##post1 treated##`x'##post1 treated##`x'##post2 treated##`x' i.year_season   ///
	   	   	   	   	    av_temp_era5_s elev_temp temp_sq  ///
	if  year<=2015 & season==`s' & elevation>1000 &  (crop_mask_v2>50|rangeland_mask_v2>50|crop_range_mask>50), fe vce(cluster mws_pixel_id_1km)

   outreg2 using "$TABLES/DID_results_`x'_`var'_2phases_season`s'.tex", ///
   keep(1.treated#1.post1 1.treated#1.post2 1.treated#1.`x'#1.post1 1.treated#1.`x'#1.post2 1.post1 1.post2  1.`x'#1.post1 1.`x'#1.post2 1.treated#1.`x' 1.`x' av_temp_era5_s  elev_temp temp_sq ) append label
    *-----------------------------------------------
    	 xtreg `var' ///
        treated##post1 treated##`x'##post1 treated##`x'##post2 treated##`x' i.year_season   ///
	   	   	   	   	    av_temp_era5_s elev_temp temp_sq  ///
						$soil ///
	if  year<=2015 & season==`s' & elevation>1000 &  (crop_mask_v2>50|rangeland_mask_v2>50|crop_range_mask>50), fe vce(cluster mws_pixel_id_1km)

   outreg2 using "$TABLES/DID_results_`x'_`var'_2phases_season`s'.tex", ///
   keep(1.treated#1.post1 1.treated#1.post2 1.treated#1.`x'#1.post1 1.treated#1.`x'#1.post2 1.post1 1.post2  1.`x'#1.post1 1.`x'#1.post2 1.treated#1.`x' 1.`x' av_temp_era5_s  elev_temp temp_sq ) append label
    *-----------------------------------------------
    	 xtreg `var'  ///
        treated##post1 treated##`x'##post1 treated##`x'##post2 treated##`x' i.year_season##i.aez_09   ///
	   	   	   	   	    av_temp_era5_s elev_temp temp_sq  ///
						$soil ///
	if  year<=2015 & season==`s' & elevation>1000 &  (crop_mask_v2>50|rangeland_mask_v2>50|crop_range_mask>50), fe vce(cluster mws_pixel_id_1km)

	    outreg2 using "$TABLES/DID_results_`x'_`var'_2phases_season`s'.tex", ///
   keep(1.treated#1.post1 1.treated#1.post2 1.treated#1.`x'#1.post1 1.treated#1.`x'#1.post2 1.post1 1.post2  1.`x'#1.post1 1.`x'#1.post2 1.treated#1.`x' 1.`x' av_temp_era5_s  elev_temp temp_sq ) append label
  *----------------------------------------------- 

 sum evi_pixel_s if treated==1 & post==0 & e(sample)==1
 local m3 = `r(mean)'
 global scale_post1_EVI_d = (1/(`m3'))*100 //% change formula
 global scale_post2_EVI_d = (1/(`m3'))*100 //% change formula
 
/* -------------------------------------------------------------> cannot identify severe drought in post1. 
 lincom 1.treated#1.`x'#1.post1+1.treated#1.post1
   mat evi_`x'= (nullmat(evi_`x')), (`r(estimate)')
   mat evi_`x'_ci = (nullmat(evi_`x'_ci)), (`r(lb)' \ `r(ub)')
*/

 lincom 1.treated#1.`x'#1.post2+1.treated#1.post2
   mat evi_`x' = (nullmat(evi_`x')), (.), (`r(estimate)')
   mat evi_`x'_ci = (nullmat(evi_`x'_ci)), ((.)\(.)), (`r(lb)' \ `r(ub)')
mat colnames evi_`x' = ATE ATE_post1 ATE_post2
mat colnames evi_`x'_ci = ATE ATE_post1 ATE_post2


 global s_`x'_post2_EVI = (1/(`m3'))*100 //% change formula
 global s_`x'_post1_EVI = (1/(`m3'))*100 //% change formula

 
  mat evi_no_drought = (nullmat(evi_no_drought)), (_b[1.treated#1.post1]), (_b[1.treated#1.post2])

  forvalues i=1/2{
 local lb_post`i' = _b[1.treated#1.post`i'] - invttail(e(df_r),0.025)*_se[1.treated#1.post`i'] 
 local ub_post`i' = _b[1.treated#1.post`i'] + invttail(e(df_r),0.025)*_se[1.treated#1.post`i']
 }
 
 mat evi_no_drought_ci =(nullmat(evi_no_drought_ci)), ((`lb_post1')\(`ub_post1')) , ((`lb_post2')\(`ub_post2')) 
 
 mat colnames evi_no_drought = ATE ATE_post1 ATE_post2
 mat colnames evi_no_drought_ci = ATE ATE_post1 ATE_post2

 *------------------------------------------------------------------------
 *-----------------------------------------------------------------------------
*---------------------------------------------------------------------------------------------------------------
    
  test _b[1.treated#1.post1]=_b[1.treated#1.post2]
  test _b[1.treated#1.post1]=0 ,notest 
  test _b[1.treated#1.post2]=0 , accum

}
}
}


 forvalues s=3/3{
 foreach var of varlist  evi_pixel_s  {
	  
foreach x in     drought_SPI_1    {
	xtreg `var' ///
       treated##`x'##post treated##`x' i.year_season ///
	if  year<=2015 & season==`s' & elevation>1000 &  (crop_mask_v2>50|rangeland_mask_v2>50|crop_range_mask>50), fe vce(cluster mws_pixel_id_1km)
	
outreg2 using "$TABLES/DID_results_`x'_`var'_season`s'.tex",keep(1.treated#1.`x'#1.post 1.treated#1.post 1.`x'#1.post 1.treated#1.`x' 1.`x' av_temp_era5_s  elev_temp temp_sq ) replace label
*-----------------------------------------------
 	xtreg `var' ///
       treated##`x'##post treated##`x' i.year_season   ///
	   	   	   	    av_temp_era5_s elev_temp temp_sq  ///
	if  year<=2015 & season==`s' & elevation>1000 &  (crop_mask_v2>50|rangeland_mask_v2>50|crop_range_mask>50), fe vce(cluster mws_pixel_id_1km)
	
	 outreg2 using "$TABLES/DID_results_`x'_`var'_season`s'.tex",keep(1.treated#1.`x'#1.post 1.treated#1.post 1.`x'#1.post 1.treated#1.`x' 1.`x' av_temp_era5_s  elev_temp temp_sq ) append label
*-----------------------------------------------
 	xtreg `var' ///
       treated##`x'##post treated##`x' i.year_season   ///
	   	   	   	    av_temp_era5_s elev_temp temp_sq  ///
					$soil ///
	if  year<=2015 & season==`s' & elevation>1000 &  (crop_mask_v2>50|rangeland_mask_v2>50|crop_range_mask>50), fe vce(cluster mws_pixel_id_1km)
	
	 outreg2 using "$TABLES/DID_results_`x'_`var'_season`s'.tex",keep(1.treated#1.`x'#1.post 1.treated#1.post 1.`x'#1.post 1.treated#1.`x' 1.`x' av_temp_era5_s  elev_temp temp_sq ) append label
*-----------------------------------------------
 	xtreg `var' ///
       treated##`x'##post treated##`x' i.year_season##i.aez_09   ///
	   	   	   	    av_temp_era5_s elev_temp temp_sq  ///
					$soil ///
	if  year<=2015 & season==`s' & elevation>1000 &  (crop_mask_v2>50|rangeland_mask_v2>50|crop_range_mask>50), fe vce(cluster mws_pixel_id_1km)
	  	 
	outreg2 using "$TABLES/DID_results_`x'_`var'_season`s'.tex",keep(1.treated#1.`x'#1.post 1.treated#1.post 1.`x'#1.post 1.treated#1.`x' 1.`x' av_temp_era5_s  elev_temp temp_sq ) append label
	
*============
	
	 sum evi_pixel_s if treated==1 & post ==0 & e(sample)==1
 local m1 = `r(mean)'
 
  global s_`x'_post_EVI = (1/(`m1'))*100 //% change formula
  global s_`x'_EVI = (1/(`m1'))*100 //% change formula
  global s_no_`x'_EVI = (1/(`m1'))*100 //% change formula

  
  lincom 1.treated#1.`x'#1.post+1.treated#1.post
   mat EVI_d_post=`r(estimate)', `r(se)', `r(p)', `r(lb)', `r(ub)'
   mat evi_`x' = `r(estimate)'
   mat evi_`x'_ci =  (`r(lb)' \ `r(ub)')

  mat evi_no_drought = (_b[1.treated#1.post])

 local lb_post = _b[1.treated#1.post] - invttail(e(df_r),0.025)*_se[1.treated#1.post] 
 local ub_post = _b[1.treated#1.post] + invttail(e(df_r),0.025)*_se[1.treated#1.post]
 
 mat evi_no_drought_ci = ((`lb_post')\(`ub_post')) 
 
  
*-----------------------------------------------	
*=======================================================

	 xtreg `var' ///
        treated##post1 treated##`x'##post1 treated##`x'##post2 treated##`x' i.year_season ///
	if  year<=2015 & season==`s' & elevation>1000 &  (crop_mask_v2>50|rangeland_mask_v2>50|crop_range_mask>50), fe vce(cluster mws_pixel_id_1km)

   outreg2 using "$TABLES/DID_results_`x'_`var'_2phases_season`s'.tex", ///
   keep(1.treated#1.post1 1.treated#1.post2 1.treated#1.`x'#1.post1 1.treated#1.`x'#1.post2 1.post1 1.post2  1.`x'#1.post1 1.`x'#1.post2 1.treated#1.`x' 1.`x' av_temp_era5_s  elev_temp temp_sq ) append label
*-----------------------------------------------   
   	 xtreg `var' ///
        treated##post1 treated##`x'##post1 treated##`x'##post2 treated##`x' i.year_season   ///
	   	   	   	   	    av_temp_era5_s elev_temp temp_sq  ///
	if  year<=2015 & season==`s' & elevation>1000 &  (crop_mask_v2>50|rangeland_mask_v2>50|crop_range_mask>50), fe vce(cluster mws_pixel_id_1km)

   outreg2 using "$TABLES/DID_results_`x'_`var'_2phases_season`s'.tex", ///
   keep(1.treated#1.post1 1.treated#1.post2 1.treated#1.`x'#1.post1 1.treated#1.`x'#1.post2 1.post1 1.post2  1.`x'#1.post1 1.`x'#1.post2 1.treated#1.`x' 1.`x' av_temp_era5_s  elev_temp temp_sq ) append label
    *-----------------------------------------------
    	 xtreg `var' ///
        treated##post1 treated##`x'##post1 treated##`x'##post2 treated##`x' i.year_season   ///
	   	   	   	   	    av_temp_era5_s elev_temp temp_sq  ///
						$soil ///
	if  year<=2015 & season==`s' & elevation>1000 &  (crop_mask_v2>50|rangeland_mask_v2>50|crop_range_mask>50), fe vce(cluster mws_pixel_id_1km)

   outreg2 using "$TABLES/DID_results_`x'_`var'_2phases_season`s'.tex", ///
   keep(1.treated#1.post1 1.treated#1.post2 1.treated#1.`x'#1.post1 1.treated#1.`x'#1.post2 1.post1 1.post2  1.`x'#1.post1 1.`x'#1.post2 1.treated#1.`x' 1.`x' av_temp_era5_s  elev_temp temp_sq ) append label
    *-----------------------------------------------
    	 xtreg `var'  ///
        treated##post1 treated##`x'##post1 treated##`x'##post2 treated##`x' i.year_season##i.aez_09   ///
	   	   	   	   	    av_temp_era5_s elev_temp temp_sq  ///
						$soil ///
	if  year<=2015 & season==`s' & elevation>1000 &  (crop_mask_v2>50|rangeland_mask_v2>50|crop_range_mask>50), fe vce(cluster mws_pixel_id_1km)

	    outreg2 using "$TABLES/DID_results_`x'_`var'_2phases_season`s'.tex", ///
   keep(1.treated#1.post1 1.treated#1.post2 1.treated#1.`x'#1.post1 1.treated#1.`x'#1.post2 1.post1 1.post2  1.`x'#1.post1 1.`x'#1.post2 1.treated#1.`x' 1.`x' av_temp_era5_s  elev_temp temp_sq ) append label
  *----------------------------------------------- 

 sum evi_pixel_s if treated==1 & post==0 & e(sample)==1
 local m3 = `r(mean)'
 global scale_post1_EVI_d = (1/(`m3'))*100 //% change formula
 global scale_post2_EVI_d = (1/(`m3'))*100 //% change formula
 

 lincom 1.treated#1.`x'#1.post1+1.treated#1.post1
   mat evi_`x'= (nullmat(evi_`x')), (`r(estimate)')
   mat evi_`x'_ci = (nullmat(evi_`x'_ci)), (`r(lb)' \ `r(ub)')



 lincom 1.treated#1.`x'#1.post2+1.treated#1.post2
   mat evi_`x' = (nullmat(evi_`x')), (`r(estimate)')
   mat evi_`x'_ci = (nullmat(evi_`x'_ci)), (`r(lb)' \ `r(ub)')
mat colnames evi_`x' = ATE ATE_post1 ATE_post2
mat colnames evi_`x'_ci = ATE ATE_post1 ATE_post2


 global s_`x'_post2_EVI = (1/(`m3'))*100 //% change formula
 global s_`x'_post1_EVI = (1/(`m3'))*100 //% change formula

 
  mat evi_no_drought = (nullmat(evi_no_drought)), (_b[1.treated#1.post1]), (_b[1.treated#1.post2])

  forvalues i=1/2{
 local lb_post`i' = _b[1.treated#1.post`i'] - invttail(e(df_r),0.025)*_se[1.treated#1.post`i'] 
 local ub_post`i' = _b[1.treated#1.post`i'] + invttail(e(df_r),0.025)*_se[1.treated#1.post`i']
 }
 
 mat evi_no_drought_ci =(nullmat(evi_no_drought_ci)), ((`lb_post1')\(`ub_post1')) , ((`lb_post2')\(`ub_post2')) 
 
 mat colnames evi_no_drought = ATE ATE_post1 ATE_post2
 mat colnames evi_no_drought_ci = ATE ATE_post1 ATE_post2

 *------------------------------------------------------------------------
 *-----------------------------------------------------------------------------
*---------------------------------------------------------------------------------------------------------------
    
  test _b[1.treated#1.post1]=_b[1.treated#1.post2]
  test _b[1.treated#1.post1]=0 ,notest 
  test _b[1.treated#1.post2]=0 , accum

}
 }
 }

 *=============================================================================================================================
 * Robustness check - Alternative land cover map *
 *----------------------------*
 
  forvalues s=3/3{
 foreach var of varlist  evi_pixel_s  {
	  	  
	   xtset mws_pixel_id_1km year_season  
			
  *###--------------------------------
  
xtreg `var' ///
       treated##post i.year_season ///
			 if  year<=2015 & season==`s' & elevation>1000 &  (crop_frac>0.50|range_frac>0.50|crop_range_frac>0.50), fe vce(cluster mws_pixel_id_1km)
			*------------------------------------------------
   outreg2 using "$TABLES/DID_results_`var'_lcover2_season`s'.tex", ///
   keep(1.treated#1.post total_precip_chirpsv20_s av_temp_era5_s total_precip_prev_s elev_temp temp_sq precip_sq elev_preci) replace label
   
xtreg `var' ///
       treated##post i.year_season total_precip_chirpsv20_s av_temp_era5_s total_precip_prev_s ///
							 elev_temp elev_preci temp_sq precip_sq  ///
			 if  year<=2015 & season==`s' & elevation>1000 &  (crop_frac>0.50|range_frac>0.50|crop_range_frac>0.50), fe vce(cluster mws_pixel_id_1km)
			*------------------------------------------------
   outreg2 using "$TABLES/DID_results_`var'_lcover2_season`s'.tex", ///
   keep(1.treated#1.post total_precip_chirpsv20_s av_temp_era5_s total_precip_prev_s elev_temp temp_sq precip_sq elev_preci) append label

xtreg `var' ///
       treated##post i.year_season total_precip_chirpsv20_s av_temp_era5_s total_precip_prev_s ///
							 elev_temp elev_preci temp_sq precip_sq  ///
							 $soil ///
			 if  year<=2015 & season==`s' & elevation>1000 &  (crop_frac>0.50|range_frac>0.50|crop_range_frac>0.50), fe vce(cluster mws_pixel_id_1km)
			*------------------------------------------------
   outreg2 using "$TABLES/DID_results_`var'_lcover2_season`s'.tex", ///
   keep(1.treated#1.post total_precip_chirpsv20_s av_temp_era5_s total_precip_prev_s elev_temp temp_sq precip_sq elev_preci) append label
	
 xtreg `var' ///
       treated##post i.year_season##i.aez_09 total_precip_chirpsv20_s av_temp_era5_s total_precip_prev_s ///
							 elev_temp elev_preci temp_sq precip_sq  ///
							 $soil ///
			 if  year<=2015 & season==`s' & elevation>1000 & (crop_frac>0.50|range_frac>0.50|crop_range_frac>0.50), fe vce(cluster mws_pixel_id_1km)
			*------------------------------------------------
   outreg2 using "$TABLES/DID_results_`var'_lcover2_season`s'.tex", ///
   keep(1.treated#1.post total_precip_chirpsv20_s av_temp_era5_s total_precip_prev_s elev_temp temp_sq precip_sq elev_preci) append label addnote(treatment group mean in post period == `m')
 

*=====================================
* Dynamic program effects *

xtreg `var' ///
       treated##post1 treated##post2 i.year_season ///
			 if  year<=2015 & season==`s' & elevation>1000 &  (crop_frac>0.50|range_frac>0.50|crop_range_frac>0.50), fe vce(cluster mws_pixel_id_1km)
			*------------------------------------------------
  outreg2 using "$TABLES/DID_results_`var'_2phases_lcover2_season`s'.tex", ///
  keep(1.treated#1.post1 1.treated#1.post2 total_precip_chirpsv20_s av_temp_era5_s total_precip_prev_s elev_temp temp_sq precip_sq elev_preci) replace label

xtreg `var' ///
       treated##post1 treated##post2 i.year_season total_precip_chirpsv20_s av_temp_era5_s total_precip_prev_s ///
							 elev_temp elev_preci temp_sq precip_sq  ///
			 if  year<=2015 & season==`s' & elevation>1000 &  (crop_frac>0.50|range_frac>0.50|crop_range_frac>0.50), fe vce(cluster mws_pixel_id_1km)
			*------------------------------------------------
  outreg2 using "$TABLES/DID_results_`var'_2phases_lcover2_season`s'.tex", ///
  keep(1.treated#1.post1 1.treated#1.post2 total_precip_chirpsv20_s av_temp_era5_s total_precip_prev_s elev_temp temp_sq precip_sq elev_preci) append label

xtreg `var' ///
       treated##post1 treated##post2 i.year_season total_precip_chirpsv20_s av_temp_era5_s total_precip_prev_s ///
							 elev_temp elev_preci temp_sq precip_sq  ///
							 $soil ///
			 if  year<=2015 & season==`s' & elevation>1000 &  (crop_frac>0.50|range_frac>0.50|crop_range_frac>0.50), fe vce(cluster mws_pixel_id_1km)
			*------------------------------------------------
  outreg2 using "$TABLES/DID_results_`var'_2phases_lcover2_season`s'.tex", ///
  keep(1.treated#1.post1 1.treated#1.post2 total_precip_chirpsv20_s av_temp_era5_s total_precip_prev_s elev_temp temp_sq precip_sq elev_preci) append label

  
xtreg `var' ///
       treated##post1 treated##post2 i.year_season##i.aez_09 total_precip_chirpsv20_s av_temp_era5_s total_precip_prev_s ///
							 elev_temp elev_preci temp_sq precip_sq  ///
							 $soil ///
			 if  year<=2015 & season==`s' & elevation>1000 &  (crop_frac>0.50|range_frac>0.50|crop_range_frac>0.50), fe vce(cluster mws_pixel_id_1km)
			*------------------------------------------------
	    outreg2 using "$TABLES/DID_results_`var'_2phases_lcover2_season`s'.tex", ///
  keep(1.treated#1.post1 1.treated#1.post2 total_precip_chirpsv20_s av_temp_era5_s total_precip_prev_s elev_temp temp_sq precip_sq elev_preci) append label 
	
 
*#################################################################################################################

 *=========================================================================
 *  DROUGHT ANALYSIS 
 *=========================================================================
 	   xtset mws_pixel_id_1km year_season 

	foreach x in  drought_SPI_2   {
	xtreg `var' ///
       treated##`x'##post treated##`x' i.year_season ///
	if  year<=2015 & season==`s' & elevation>1000 &  (crop_frac>0.50|range_frac>0.50|crop_range_frac>0.50), fe vce(cluster mws_pixel_id_1km)
	
	 outreg2 using "$TABLES/DID_results_`x'_`var'_lcover2_season`s'.tex",keep(1.treated#1.`x'#1.post 1.treated#1.post 1.`x'#1.post 1.treated#1.`x' 1.`x' av_temp_era5_s  elev_temp temp_sq ) replace label
*-----------------------------------------------
 	xtreg `var' ///
       treated##`x'##post treated##`x' i.year_season   ///
	   	   	   	    av_temp_era5_s elev_temp temp_sq  ///
	if  year<=2015 & season==`s' & elevation>1000 &  (crop_frac>0.50|range_frac>0.50|crop_range_frac>0.50), fe vce(cluster mws_pixel_id_1km)
	
	 outreg2 using "$TABLES/DID_results_`x'_`var'_lcover2_season`s'.tex",keep(1.treated#1.`x'#1.post 1.treated#1.post 1.`x'#1.post 1.treated#1.`x' 1.`x' av_temp_era5_s  elev_temp temp_sq ) append label
*-----------------------------------------------
 	xtreg `var' ///
       treated##`x'##post treated##`x' i.year_season   ///
	   	   	   	    av_temp_era5_s elev_temp temp_sq  ///
					$soil ///
	if  year<=2015 & season==`s' & elevation>1000 &  (crop_frac>0.50|range_frac>0.50|crop_range_frac>0.50), fe vce(cluster mws_pixel_id_1km)
	
	 outreg2 using "$TABLES/DID_results_`x'_`var'_lcover2_season`s'.tex",keep(1.treated#1.`x'#1.post 1.treated#1.post 1.`x'#1.post 1.treated#1.`x' 1.`x' av_temp_era5_s  elev_temp temp_sq ) append label
*-----------------------------------------------
 	xtreg `var' ///
       treated##`x'##post treated##`x' i.year_season##i.aez_09   ///
	   	   	   	    av_temp_era5_s elev_temp temp_sq  ///
					$soil ///
	if  year<=2015 & season==`s' & elevation>1000 &  (crop_frac>0.50|range_frac>0.50|crop_range_frac>0.50), fe vce(cluster mws_pixel_id_1km)
	  	 
	outreg2 using "$TABLES/DID_results_`x'_`var'_lcover2_season`s'.tex",keep(1.treated#1.`x'#1.post 1.treated#1.post 1.`x'#1.post 1.treated#1.`x' 1.`x' av_temp_era5_s  elev_temp temp_sq ) append label
	
*-----------------------------------------------	
*=======================================================

	 xtreg `var' ///
        treated##post1 treated##`x'##post1 treated##`x'##post2 treated##`x' i.year_season ///
	if  year<=2015 & season==`s' & elevation>1000 &  (crop_frac>0.50|range_frac>0.50|crop_range_frac>0.50), fe vce(cluster mws_pixel_id_1km)

   outreg2 using "$TABLES/DID_results_`x'_`var'_2phases_lcover2_season`s'.tex", ///
   keep(1.treated#1.post1 1.treated#1.post2 1.treated#1.`x'#1.post1 1.treated#1.`x'#1.post2 1.post1 1.post2  1.`x'#1.post1 1.`x'#1.post2 1.treated#1.`x' 1.`x' av_temp_era5_s  elev_temp temp_sq ) append label
*-----------------------------------------------   
   	 xtreg `var' ///
        treated##post1 treated##`x'##post1 treated##`x'##post2 treated##`x' i.year_season   ///
	   	   	   	   	    av_temp_era5_s elev_temp temp_sq  ///
	if  year<=2015 & season==`s' & elevation>1000 &  (crop_frac>0.50|range_frac>0.50|crop_range_frac>0.50), fe vce(cluster mws_pixel_id_1km)

   outreg2 using "$TABLES/DID_results_`x'_`var'_2phases_lcover2_season`s'.tex", ///
   keep(1.treated#1.post1 1.treated#1.post2 1.treated#1.`x'#1.post1 1.treated#1.`x'#1.post2 1.post1 1.post2  1.`x'#1.post1 1.`x'#1.post2 1.treated#1.`x' 1.`x' av_temp_era5_s  elev_temp temp_sq ) append label
    *-----------------------------------------------
    	 xtreg `var' ///
        treated##post1 treated##`x'##post1 treated##`x'##post2 treated##`x' i.year_season   ///
	   	   	   	   	    av_temp_era5_s elev_temp temp_sq  ///
						$soil ///
	if  year<=2015 & season==`s' & elevation>1000 &  (crop_frac>0.50|range_frac>0.50|crop_range_frac>0.50), fe vce(cluster mws_pixel_id_1km)

   outreg2 using "$TABLES/DID_results_`x'_`var'_2phases_lcover2_season`s'.tex", ///
   keep(1.treated#1.post1 1.treated#1.post2 1.treated#1.`x'#1.post1 1.treated#1.`x'#1.post2 1.post1 1.post2  1.`x'#1.post1 1.`x'#1.post2 1.treated#1.`x' 1.`x' av_temp_era5_s  elev_temp temp_sq ) append label
    *-----------------------------------------------
    	 xtreg `var'  ///
        treated##post1 treated##`x'##post1 treated##`x'##post2 treated##`x' i.year_season##i.aez_09   ///
	   	   	   	   	    av_temp_era5_s elev_temp temp_sq  ///
						$soil ///
	if  year<=2015 & season==`s' & elevation>1000 &  (crop_frac>0.50|range_frac>0.50|crop_range_frac>0.50), fe vce(cluster mws_pixel_id_1km)

	    outreg2 using "$TABLES/DID_results_`x'_`var'_2phases_lcover2_season`s'.tex", ///
   keep(1.treated#1.post1 1.treated#1.post2 1.treated#1.`x'#1.post1 1.treated#1.`x'#1.post2 1.post1 1.post2  1.`x'#1.post1 1.`x'#1.post2 1.treated#1.`x' 1.`x' av_temp_era5_s  elev_temp temp_sq ) append label
  *----------------------------------------------- 

}
}
}

 
*#################################################################################################################
* ############################################################# SIF, GOSIF-GPP, 5-km EVI ############################################

	use "$INPUT/SLMP_SIF_pixel_level", replace 

gen treated=(donor=="WBI")
lab define treated 1 "Treated[SLMP1]" 0 "Control[SLMP2]"
la val treated treated 

 tostring year, gen(year_s)
 tostring month, gen(month_s)
 
 
* Creating season variables to average satelite indicators over - following FAO's seasonality classification for Ethiopia
 gen season=.
	replace season=1 if (month>=10 | month==1)
	replace season=2 if  month>=2 & month<6
	replace season=3 if  month>=6 & month<10
lab define season 1 "Dry Season" 2 "Pre-rainy season" 3 "Rainy season"

la val season season 
  tostring season, gen(season_s)
  gen year_season_s=year_s+season_s
  destring year_season_s, gen(year_season)

  local vars  total_precip_chirpsv20 av_temp_era5
  foreach var of varlist `vars' {
  bys mws_pixel_id_5km year_season: egen `var'_s=mean(`var') // averaging at the year_season level
  }
  
  rename gosif_gpp_ gosif_gpp  
  
 foreach var of varlist sif gosif_gpp evi_5km{
 
preserve 
		 sum `var'
		   gen nm=cond(`var'!=.,0,1)
		   bys mws_pixel_id_5km year_season: egen missing=total(nm)
		   drop if missing!=0
		   
		  bys mws_pixel_id_5km year_season: egen `var'_s=mean(`var') // averaging at the year_season level
		  bys mws_pixel_id_5km year_season: egen `var'_max_s=max(`var') // max of the the year_season level

		  keep mws_pixel_id_5km year_season `var'_s  `var'_max_s  
		  bys mws_pixel_id_5km year_season: keep if _n==1
		  tempfile average_`var'
		  save `average_`var'', replace 
 restore 
  
    merge m:1 mws_pixel_id_5km year_season using `average_`var'', nogen
	}

	  gen crop_range_mask=(crop_mask_v2+rangeland_mask_v2)
	  gen crop_range_frac=(crop_frac+range_frac)

*============================================================
        *# GENERATING DROUGHT VARIABLE #* 
   
   gen drought_SPI_1=(SPI<=-1) if SPI!=.
   gen drought_SPI_15=(SPI<=-1.5) if SPI!=.
   gen drought_SPI_2=(SPI<=-2) if SPI!=.

     foreach var of varlist drought* {
   bys mws_pixel_id_5km year_season: egen max_`var'=max(`var')
   replace `var'= max_`var'
   drop max_`var'
  }
  
  replace sif_s=. if year==2002 //no complete series for the rainy season in this year so eliminating it
  
  
  
      * ==== GRAPH OF DROUGHT OVER TIME ==== *

	   preserve
		keep if elevation>1000 & (crop_mask_v2>50|rangeland_mask_v2>50|crop_range_mask>50)
		bys mws_pixel_id_5km year_season: keep if _n==1
		
		foreach var of varlist drought*{
			bys year_season: egen prop_`var'_1 = mean(`var') if treated==1
			bys year_season: egen prop_`var'_0 = mean(`var') if treated==0
		
		
 twoway (line  prop_`var'_1  year if treated==1 & season==3 & year<=2015, lcolor(orange red) lwidth(medthick) ///
         ylabel(0(0.2)1, format(%03.2f) labsize(medium) glpattern(dot) glcolor(gray%90) glwidth(thick) ) ///
		 legend(lab(1 "SLMP1") lab(2 "SLMP2") col(2) position(6) size(medium)) ///
		 xlabel(#16,  glpattern(dot) glcolor(gray%90) glwidth(thick) labsize(medium)) xtitle("", height(10) ) ///
		 ytitle("Prop. of pixels that suffered a drought event", height(10) size(medlarge)) ) ///
		 (line  prop_`var'_0 year if treated==0 & season==3 & year<=2015, lcolor(black%50) lpattern(dash) lwidth(medthick))
  graph export "$GRAPHS/SI_Fig1_lpolyreg_`var'_5km_season3.tif", replace
		}
     restore 
	
 *============================================================================================
 *##### SUMMARY STATISTICS ######*
 
* Time variant vars *

 preserve
  keep if season==3 & evi_5km_s!=. //just pixels that form part of the estimation model
  bys mws_pixel_id_5km year_season: keep if _n==1
  keep if year<=2015 & year>=2002
 
	*---------------------------------
	
* Precipitation and temperature graphs *

 rename total_precip_chirpsv20_s t_precip
 local vars t_precip av_temp_era5_s 
 levelsof year_season, local(year_season)

 forvalues i=0/1{
   foreach var of varlist `vars'{
		gen lb_`var'_`i' = .
		gen ub_`var'_`i' = .
		
foreach l of local year_season {
    ci mean `var' if year_season == `l' & treated==`i'  & elevation>1000 & (crop_mask_v2>50|rangeland_mask_v2>50|crop_range_mask>50)
    replace lb_`var'_`i' = `r(lb)' if year_season == `l' & treated==`i'
    replace ub_`var'_`i' = `r(ub)' if year_season == `l' & treated==`i'
}
}
}

  foreach var of varlist `vars'{
 bys year_season: egen  mean_`var'_1=mean(`var') if treated==1  & elevation>1000 & (crop_mask_v2>50|rangeland_mask_v2>50|crop_range_mask>50)
 bys year_season: egen  mean_`var'_0=mean(`var') if treated==0  & elevation>1000 & (crop_mask_v2>50|rangeland_mask_v2>50|crop_range_mask>50)

}

 twoway (rarea lb_t_precip_1 ub_t_precip_1 year, alcolor(black%20) fcolor(orange red) fintensity(inten10) msize(*2)) ///
        (line  mean_t_precip_1 year if treated==1,  lpattern(solid) lcolor(orange red) lwidth(medthick) ///
		ylabel(150(25)250, glpattern(dot) glcolor(gray%90) glwidth(thick)  labsize(medium)) ///
		legend(lab(1 "CI (95%)") lab(2 "SLMP1") lab(3 "SLMP2") lab(4 "CI (95%)") col(4) position(6) size(medsmall)) ///
		xlabel(#16, labsize(medium) glpattern(dot) glcolor(gray%90) glwidth(thick) ) ///
		xtitle("", height(10)) ytitle("Seasonal average monthly precipitation (mm)", size(medlarge) height(10)) ) ///
        (rarea lb_t_precip_0 ub_t_precip_0 year, alcolor(black%20) fcolor(white)  fintensity(inten10) msize(*2)) ///
		(line  mean_t_precip_0 year if treated==0,  lcolor(black) lpattern(dash) lwidth(medthick) ) 
 
    graph export "$GRAPHS/SI_Figure7_leftpanel_Precipitation_over_time.tif", replace

 twoway (rarea lb_av_temp_era5_s_1 ub_av_temp_era5_s_1 year, alcolor(black%20) fcolor(orange red) fintensity(inten30)  msize(*2)) ///
        (line  mean_av_temp_era5_s_1 year if treated==1,  lpattern(solid) lcolor(orange red) lwidth(medthick) ///
		ylabel(, glpattern(dot) glcolor(gray%90) glwidth(thick)  labsize(medium)) ///
		legend(lab(1 "CI (95%)") lab(2 "SLMP1") lab(3 "SLMP2") lab(4 "CI (95%)") col(4) position(6) size(medsmall)) ///
		xlabel(#16, labsize(medium) glpattern(dot) glcolor(gray%90) glwidth(thick) ) ///
		xtitle("", height(10)) ytitle("Seasonal average monthly temperature ( {sup:{c 176}}C)", size(medlarge) height(10)) ) ///
        (rarea lb_av_temp_era5_s_0 ub_av_temp_era5_s_0 year, alcolor(black%20) fcolor(white%30)  fintensity(inten30) msize(*2)) ///
		(line  mean_av_temp_era5_s_0 year if treated==0,  lcolor(black) lpattern(dash) lwidth(medthick) ) 
 
    graph export "$GRAPHS/SI_Figure7_rightpanel_Temperature_over_time.tif", replace

restore 

   *============================================================================================
***#############################
	 *---------------------------------------------
 * 1. Cheking for paralell trends assumption 
 *--------------------------------------------
		*1.1. Visual approach - local polynomial non parametric regression of EVI/SIF on time by treatment group
  
  *----------
  *---------- 

 preserve
  keep if season==3 & year<=2015 & year>=2002
    bys year_season: gen time_cons=1 if _n==1 
  replace time_cons=sum(time_cons)
  replace time_cons=time_cons-1

  
    twoway (lpolyci  gosif_gpp_s time_cons if treated==1 & elevation>1000 & (crop_mask_v2>50|rangeland_mask_v2>50|crop_range_mask>50) ,   bwidth(0.5) ///
	legend(order(2 "SLMP1" 4 "SLMP2") col(2) position(6) height(5) size(medium)) xline(8.5 13) ///
	ytitle("Seasonal average GOSIF-GPP (g C m{sup:-2} month{sup:-1})", height(8) size(medlarge)) lwidth(medthick) fcolor(orange red) fintensity(inten10) ///
	xlabel(0 "2002" 1 "2003" 2 "2004" 3 "2005" 4 "2006" 5 "2007" 6 "2008" 7 "2009" 8 "2010" 9 "2011" 10 "2012" 11 "2013" 12 "2014" 13 "2015",    ///
	 glpattern(dot) glcolor(gray%90) glwidth(thick)  labsize(medium)) ylabel(, glpattern(dot) glcolor(gray%90) glwidth(thick)  labsize(medium)))   ///
	(lpolyci  gosif_gpp_s time_cons if treated==0 & elevation>1000 & (crop_mask_v2>50|rangeland_mask_v2>50|crop_range_mask>50), bwidth(0.5) ///
	  lcolor(black)  alcolor(black%20) ciplot(rline) alpattern(solid) fcolor(orange red) fintensity(inten10) lpattern(shortdash) xtitle("") lwidth(medthick))
  graph export "$GRAPHS/Fig2_rightpanel_GOSIF_GPP_lpolyreg_season3.tif", replace
  
  
   twoway (lpolyci  sif_s time_cons if treated==1 & elevation>1000 & (crop_mask_v2>50|rangeland_mask_v2>50|crop_range_mask>50) , bwidth(0.5) ///
	legend(order(2 "SLMP1" 4 "SLMP2") col(2) position(6) height(5) size(medium)) xline(8.5 13) ///
	ytitle("SIF (W m{sup:-2} {&mu}m{sup:-1} sr{sup:-1})", height(8) size(medlarge)) lwidth(medthick) fcolor(orange red) fintensity(inten10) ///
	xlabel(0 "2002" 1 "2003" 2 "2004" 3 "2005" 4 "2006" 5 "2007" 6 "2008" 7 "2009" 8 "2010" 9 "2011" 10 "2012" 11 "2013" 12 "2014" 13 "2015",    ///
	  glpattern(dot) glcolor(gray%90) glwidth(thick)  labsize(medium)) ylabel(0.70(0.05)1, glpattern(dot) glcolor(gray%90) glwidth(thick) format(%03.2f) labsize(medium)))   ///
	(lpolyci  sif_s time_cons if treated==0 & elevation>1000 & (crop_mask_v2>50|rangeland_mask_v2>50|crop_range_mask>50),  bwidth(0.5) ///
	  lcolor(black)  alcolor(black%20) ciplot(rline) alpattern(solid) fcolor(orange red) fintensity(inten10) lpattern(shortdash) xtitle("") lwidth(medthick))
  graph export "$GRAPHS/SI_Fig2_rightpanel_SIF_lpolyreg_season3.tif", replace

   twoway (lpolyci  evi_5km_s time_cons if treated==1 & elevation>1000 & (crop_mask_v2>50|rangeland_mask_v2>50|crop_range_mask>50) ,  bwidth(0.5) ///
	legend(order(2 "SLMP1" 4 "SLMP2") col(2) position(6) height(5) size(medium)) xline(8.5 13) ///
	ytitle("EVI", height(8) size(medlarge)) lwidth(medthick) fcolor(orange red) fintensity(inten10) ///
	xlabel(0 "2002" 1 "2003" 2 "2004" 3 "2005" 4 "2006" 5 "2007" 6 "2008" 7 "2009" 8 "2010" 9 "2011" 10 "2012" 11 "2013" 12 "2014" 13 "2015",  ///
	 glpattern(dot) glcolor(gray%90) glwidth(thick)  labsize(medium)) ylabel(0.32(0.02)0.46, glpattern(dot) glcolor(gray%90) glwidth(thick)  format(%03.2f) labsize(medium)))   ///
	(lpolyci  evi_5km_s time_cons if treated==0 & elevation>1000 & (crop_mask_v2>50|rangeland_mask_v2>50|crop_range_mask>50), bwidth(0.5) ///
	  lcolor(black)  alcolor(black%20) ciplot(rline) alpattern(solid) fcolor(orange red) fintensity(inten10) lpattern(shortdash) xtitle("") lwidth(medthick))
  graph export "$GRAPHS/SI_Fig2_leftpanel_EVI_5km_lpolyreg_season3.tif", replace
  
  restore 
  
 *-------------------------------------------------
 
    gen post=(year>2010) // based on budget disbursement 
	 gen post1=(year>2010 & year<=2013)
	 gen post2=(year>=2014 & year<=2015)

   
 *=======================================================
 *=========================================================================
 *  SPECIFICATION for parallel trends check : EVENT STUDY
 *=========================================================================
  
   bys mws_pixel_id_5km year_season: keep if _n==1
   drop if mws_pixel_id_5km==.

	  gen elev_preci=total_precip_chirpsv20_s*elevation
		 gen elev_temp=av_temp_era5_s*elevation
		  gen temp_sq = av_temp_era5_s^2
		  gen precip_sq = total_precip_chirpsv20_s^2
		  gen precip_wate= total_precip_chirpsv20_s*waterholdi

		  
  bysort mws_pixel_id_5km (year_season) : gen total_precip_prev_s = (total_precip_chirpsv20_s[_n-1])

  	  foreach var of varlist af* {
	  gen prec_`var'=`var'*total_precip_chirpsv20_s
	  gen temp_`var'=`var'*av_temp_era5_s	  
	  }

 global soil "prec_af250m_nut prec_af250m_n_1 prec_af250m_n_2  prec_af_pH515cm prec_af_ORCDRC1 prec_af_EC515cm  prec_af_Depthto prec_af_CEC515c precip_wate"

*------------------

tab year_season if year<=2015 & season==3, gen(period)
 foreach var of varlist period*{
  replace `var'=`var'*treated 
 }
 
 la var period1 "period -9 (year 2002)"
 la var period2 "period -8  (year 2003)"
 la var period3 "period -7  (year 2004)"
 la var period4 "period -6  (year 2005)"
 la var period5 "period -5  (year 2006)"
 la var period6 "period -4  (year 2007)"
 la var period7 "period -3  (year 2008)"
 la var period8 "period -2  (year 2009)"
 la var period10 "period +1  (year 2011)"
 la var period11 "period +2  (year 2012)"
 la var period12 "period +3  (year 2013)"
 la var period13 "period +4  (year 2014)"
 la var period14 "period +5  (year 2015)"

  constraint define 1  period1 + period2 + period3 + period4 + period5 ///
 + period6 + period7 + period8   = 0
	
 foreach var of varlist evi_5km_s  sif_s gosif_gpp_s    {

	cnsreg `var' ///
		 period1-period10 period11 period12-period14 ///
				 i.mws_pixel_id_5km i.year_season ///
			 if  year<=2015  & season==3 & elevation>1000 & (crop_mask_v2>50|rangeland_mask_v2>50|crop_range_mask>50), noconstant constraint(1) vce(cluster mws_pixel_id_5km) 
	outreg2 using "$TABLES/Event_Study_results_`var'_Millermethod.tex", keep(period*) replace label
		 
		predict fitted if e(sample)==1
		corr `var' fitted if e(sample)==1
		drop fitted 
		
*----------------------------------------------------
  	cnsreg `var' ///
		 period1-period10 period11 period12-period14 ///
               total_precip_chirpsv20_s av_temp_era5_s total_precip_prev_s ///
				elev_temp elev_preci temp_sq precip_sq precip_wate  /// 
				 i.mws_pixel_id_5km i.year_season ///
			 if  year<=2015  & season==3 & elevation>1000 & (crop_mask_v2>50|rangeland_mask_v2>50|crop_range_mask>50), noconstant constraint(1) vce(cluster mws_pixel_id_5km) 
    outreg2 using "$TABLES/Event_Study_results_`var'_Millermethod.tex", keep(period*) append label
   
   	 	predict fitted if e(sample)==1
		corr `var' fitted if e(sample)==1
		drop fitted 
*----------------------------------------------------
  	cnsreg `var' ///
		 period1-period10 period11 period12-period14 ///
               total_precip_chirpsv20_s av_temp_era5_s total_precip_prev_s ///
				elev_temp elev_preci temp_sq precip_sq precip_wate  /// 
				$soil ///
				 i.mws_pixel_id_5km i.year_season ///
			 if  year<=2015  & season==3 & elevation>1000 & (crop_mask_v2>50|rangeland_mask_v2>50|crop_range_mask>50), noconstant constraint(1) vce(cluster mws_pixel_id_5km) 
   	 outreg2 using "$TABLES/Event_Study_results_`var'_Millermethod.tex", keep(period*) append label
	 
	 	predict fitted if e(sample)==1
		corr `var' fitted if e(sample)==1
		drop fitted 
*----------------------------------------------------
  	cnsreg `var' ///
		 period1-period10 period11 period12-period14 ///
               total_precip_chirpsv20_s av_temp_era5_s total_precip_prev_s ///
				elev_temp elev_preci temp_sq precip_sq precip_wate  /// 
				$soil ///
				 i.mws_pixel_id_5km i.year_season##i.aez_09 ///
			 if  year<=2015  & season==3 & elevation>1000 & (crop_mask_v2>50|rangeland_mask_v2>50|crop_range_mask>50), noconstant constraint(1) vce(cluster mws_pixel_id_5km) 
   	  outreg2 using "$TABLES/Event_Study_results_`var'_Millermethod.tex", keep(period*) append label
	
		predict fitted if e(sample)==1
		corr `var' fitted if e(sample)==1
		drop fitted 
*----------------------------------------------------
*----------------------------------------------------	  
   foreach i in 1 2 3 4 5 6 7 8 9  10 11  12 13 14 {
    gen beta_`var'_`i'=_b[period`i'] if period`i'==1
	gen se_`var'_`i'=_se[period`i'] if period`i'==1
   }
   
   egen beta_`var'=rowmax(beta_`var'_*)
   egen se_`var'=rowmax(se_`var'_*)	 

 }
   	
	   sort year_season

     serrbar beta_gosif_gpp_s se_gosif_gpp_s year if year<=2015, scale (1.96)  xlabel(#15) ///
   ytitle("GOSIF-GPP (g C m{sup:-2} month{sup:-1})", size(medlarge) height(8)) xtitle("") ///
   yline(0, lwidth(medthick)) mvopts(c(i)  m(D) msize(medium) mlw(medlarge) mlc(black%20) mc(gray%50)) ///
   xlabel(,  glpattern(dot) glcolor(gray%90) glwidth(thick) labsize(medium)) ylabel(, glpattern(dot) glcolor(gray%90) glwidth(thick)  labsize(medium)) ///
    lwidth(medthick) lcolor(orange_red) xline(2010.50, lcolor(orange_red) lwidth(medthick) lpattern(longdash)) 
   graph export "$GRAPHS/Fig3_panelB_gosif_gpp_Event_Study_Miller_method.tif", replace
   
    
    serrbar beta_evi_5km_s se_evi_5km_s year if year<=2015, scale (1.96)  xlabel(#15) ///
   ytitle("EVI", size(medlarge) height(8)) xtitle("") ///
   yline(0, lwidth(medthick)) mvopts(c(i)  m(D) msize(medium) mlw(medlarge) mlc(black%20) mc(gray%50)) ///
   xlabel(, glpattern(dot) glcolor(gray%90) glwidth(thick) labsize(medium)) ylabel(-0.03(0.01)0.03, glpattern(dot) glcolor(gray%90) glwidth(thick) format(%03.2f) labsize(medium)) ///
   lwidth(medthick) lcolor(orange_red) xline(2010.50, lcolor(orange_red) lwidth(medthick) lpattern(longdash)) 
   graph export "$GRAPHS/SI_Fig3_panelA_evi_5km_Event_Study_Miller_method.tif", replace
 
   	
   	serrbar beta_sif_s se_sif_s year if year<=2015, scale (1.96)  xlabel(#15) ///
   ytitle("SIF (W m{sup:-2} {&mu}m{sup:-1} sr{sup:-1})", size(medlarge) height(8)) xtitle("") ///
   yline(0, lwidth(medthick)) mvopts(c(i)  m(D) msize(medium) mlw(medlarge) mlc(black%20) mc(gray%50)) ///
   xlabel(, glpattern(dot) glcolor(gray%90) glwidth(thick) labsize(medium)) ylabel(-0.03(0.01)0.03, glpattern(dot) glcolor(gray%90) glwidth(thick) format(%03.2f) labsize(medium)) ///
   lwidth(medthick) lcolor(orange_red) xline(2010.50, lcolor(orange_red) lwidth(medthick) lpattern(longdash)) 
   graph export "$GRAPHS/SI_Fig3_panelB_sif_Event_Study_Miller_method.tif", replace
   
	drop beta_* se_*	 

*----------------------------------------------------------------------------------------------	
* Robustness check - alternative land cover map

foreach var of varlist  evi_5km_s sif_s  gosif_gpp_s    {

	cnsreg `var' ///
		 period1-period10 period11 period12-period14 ///
				 i.mws_pixel_id_5km i.year_season ///
			 if  year<=2015  & season==3 & elevation>1000 & (crop_frac>0.50|range_frac>0.50|crop_range_frac>0.50), noconstant constraint(1) vce(cluster mws_pixel_id_5km) 
	outreg2 using "$TABLES/Event_Study_results_`var'_altcover.tex", keep(period*) replace label
		 
		predict fitted if e(sample)==1
		corr `var' fitted if e(sample)==1
		drop fitted 
		
*----------------------------------------------------
  	cnsreg `var' ///
		 period1-period10 period11 period12-period14 ///
               total_precip_chirpsv20_s av_temp_era5_s total_precip_prev_s ///
				elev_temp elev_preci temp_sq precip_sq precip_wate  /// 
				 i.mws_pixel_id_5km i.year_season ///
			 if  year<=2015  & season==3 & elevation>1000 & (crop_frac>0.50|range_frac>0.50|crop_range_frac>0.50), noconstant constraint(1) vce(cluster mws_pixel_id_5km) 
    outreg2 using "$TABLES/Event_Study_results_`var'_altcover.tex", keep(period*) append label
   
   	 	predict fitted if e(sample)==1
		corr `var' fitted if e(sample)==1
		drop fitted 
*----------------------------------------------------
  	cnsreg `var' ///
		 period1-period10 period11 period12-period14 ///
               total_precip_chirpsv20_s av_temp_era5_s total_precip_prev_s ///
				elev_temp elev_preci temp_sq precip_sq precip_wate  /// 
				$soil ///
				 i.mws_pixel_id_5km i.year_season ///
			 if  year<=2015  & season==3 & elevation>1000 & (crop_frac>0.50|range_frac>0.50|crop_range_frac>0.50), noconstant constraint(1) vce(cluster mws_pixel_id_5km) 
   	 outreg2 using "$TABLES/Event_Study_results_`var'_altcover.tex", keep(period*) append label
	 
	 	predict fitted if e(sample)==1
		corr `var' fitted if e(sample)==1
		drop fitted 
*----------------------------------------------------
  	cnsreg `var' ///
		 period1-period10 period11 period12-period14 ///
               total_precip_chirpsv20_s av_temp_era5_s total_precip_prev_s ///
				elev_temp elev_preci temp_sq precip_sq precip_wate  /// 
				$soil ///
				 i.mws_pixel_id_5km i.year_season##i.aez_09 ///
			 if  year<=2015  & season==3 & elevation>1000 & (crop_frac>0.50|range_frac>0.50|crop_range_frac>0.50), noconstant constraint(1) vce(cluster mws_pixel_id_5km) 
   	  outreg2 using "$TABLES/Event_Study_results_`var'_altcover.tex", keep(period*) append label
	
		predict fitted if e(sample)==1
		corr `var' fitted if e(sample)==1
		drop fitted 
*----------------------------------------------------
*----------------------------------------------------	  
   foreach i in 1 2 3 4 5 6 7 8 9  10 11  12 13 14 {
    gen beta_`var'_`i'=_b[period`i'] if period`i'==1
	gen se_`var'_`i'=_se[period`i'] if period`i'==1
   }
   
   egen beta_`var'=rowmax(beta_`var'_*)
   egen se_`var'=rowmax(se_`var'_*)
	
 }
   	
	   sort year_season

     serrbar beta_gosif_gpp_s se_gosif_gpp_s year if year<=2015, scale (1.96)  xlabel(#15) ///
   ytitle("GOSIF-GPP (g C m{sup:-2} month{sup:-1})", size(medlarge) height(8)) xtitle("") ///
   yline(0, lwidth(medthick)) mvopts(c(i)  m(D) msize(medium) mlw(medlarge) mlc(black%20) mc(gray%50)) ///
   xlabel(,  glpattern(dot) glcolor(gray%90) glwidth(thick) labsize(medium)) ylabel(, glpattern(dot) glcolor(gray%90) glwidth(thick)  labsize(medium)) ///
    lwidth(medthick) lcolor(orange_red) xline(2010.50, lcolor(orange_red) lwidth(medthick) lpattern(longdash)) 
   graph export "$GRAPHS/SI_Fig4_panelB_gosif_gpp_Event_Study_altcover.tif", replace
   
    
    serrbar beta_evi_5km_s se_evi_5km_s year if year<=2015, scale (1.96)  xlabel(#15) ///
   ytitle("EVI", size(medlarge) height(8)) xtitle("") ///
   yline(0, lwidth(medthick)) mvopts(c(i)  m(D) msize(medium) mlw(medlarge) mlc(black%20) mc(gray%50)) ///
   xlabel(, glpattern(dot) glcolor(gray%90) glwidth(thick) labsize(medium)) ylabel(-0.03(0.01)0.03, glpattern(dot) glcolor(gray%90) glwidth(thick) format(%03.2f) labsize(medium)) ///
   lwidth(medthick) lcolor(orange_red) xline(2010.50, lcolor(orange_red) lwidth(medthick) lpattern(longdash)) 
   graph export "$GRAPHS/SI_Fig4_panelC_evi_5km_Event_Study_altcover.tif", replace
 
   	
   	serrbar beta_sif_s se_sif_s year if year<=2015, scale (1.96)  xlabel(#15) ///
   ytitle("SIF (W m{sup:-2} {&mu}m{sup:-1} sr{sup:-1})", size(medlarge) height(8)) xtitle("") ///
   yline(0, lwidth(medthick)) mvopts(c(i)  m(D) msize(medium) mlw(medlarge) mlc(black%20) mc(gray%50)) ///
   xlabel(, glpattern(dot) glcolor(gray%90) glwidth(thick) labsize(medium)) ylabel(-0.03(0.01)0.03, glpattern(dot) glcolor(gray%90) glwidth(thick) format(%03.2f) labsize(medium)) ///
   lwidth(medthick) lcolor(orange_red) xline(2010.50, lcolor(orange_red) lwidth(medthick) lpattern(longdash)) 
   graph export "$GRAPHS/SI_Fig4_panelD_sif_Event_Study_altcover.tif", replace
	
   
	drop beta_* se_*	 
	
*=====================================================

*####################################################
			
 *=========================================================================
 *  DIFFERENCE IN DIFFERENCE
 *=========================================================================
 rename drought_SPI_15 d_SPI_15
 
 forvalues s =3/3{
 foreach var of varlist  evi_5km_s sif_s gosif_gpp_s {
   
	   xtset mws_pixel_id_5km year_season 
  
  *###--------------------------------
	
xtreg `var' ///
      treated##post  i.year_season ///
			 if  year<=2015 & season==`s' & elevation>1000 & (crop_mask_v2>50|rangeland_mask_v2>50|crop_range_mask>50), fe vce(cluster mws_pixel_id_5km)
			*------------------------------------------------
	
	outreg2 using "$TABLES/DID_results_`var'_season`s'.tex", ///
	 keep(1.treated#1.post total_precip_chirpsv20_s av_temp_era5_s total_precip_prev_s elev_temp temp_sq precip_sq elev_preci) replace label
*----------------------------------------------------	 
xtreg `var' ///
      treated##post  i.year_season   total_precip_chirpsv20_s  av_temp_era5_s total_precip_prev_s ///
							elev_preci elev_temp temp_sq precip_sq precip_wate   ///
			 if  year<=2015 & season==`s' & elevation>1000 & (crop_mask_v2>50|rangeland_mask_v2>50|crop_range_mask>50), fe vce(cluster mws_pixel_id_5km)
			*------------------------------------------------
	
	outreg2 using "$TABLES/DID_results_`var'_season`s'.tex", ///
	 keep(1.treated#1.post total_precip_chirpsv20_s av_temp_era5_s total_precip_prev_s elev_temp temp_sq precip_sq elev_preci) append label
*----------------------------------------------------	
xtreg `var' ///
      treated##post  i.year_season   total_precip_chirpsv20_s  av_temp_era5_s total_precip_prev_s ///
							elev_preci elev_temp temp_sq precip_sq precip_wate   ///
							$soil ///
			 if  year<=2015 & season==`s' & elevation>1000 & (crop_mask_v2>50|rangeland_mask_v2>50|crop_range_mask>50), fe vce(cluster mws_pixel_id_5km)
			*------------------------------------------------
	
	outreg2 using "$TABLES/DID_results_`var'_season`s'.tex", ///
	 keep(1.treated#1.post total_precip_chirpsv20_s av_temp_era5_s total_precip_prev_s elev_temp temp_sq precip_sq elev_preci) append label
*----------------------------------------------------	
xtreg `var' ///
      treated##post  i.year_season##i.aez_09   total_precip_chirpsv20_s  av_temp_era5_s total_precip_prev_s ///
							elev_preci elev_temp temp_sq precip_sq precip_wate   ///
							$soil ///
			 if  year<=2015 & season==`s' & elevation>1000 & (crop_mask_v2>50|rangeland_mask_v2>50|crop_range_mask>50), fe vce(cluster mws_pixel_id_5km)
			*------------------------------------------------
  	outreg2 using "$TABLES/DID_results_`var'_season`s'.tex", ///
	 keep(1.treated#1.post total_precip_chirpsv20_s av_temp_era5_s total_precip_prev_s elev_temp temp_sq precip_sq elev_preci) append label ///
	 addnote(treatment group mean in the post period = `m')
	 
 estimates store DID_`var'

 * Coefficients for percentage change * 
 
 
 sum `var' if treated==1 & post ==0 & e(sample)==1
 local m = `r(mean)'
 
 global scale1_`var'_`s' = (1/(`m'))*100 //% change formula
 

 mat `var'_general_season`s' = _b[1.treated#1.post]


 local lb_post = _b[1.treated#1.post] - invttail(e(df_r),0.025)*_se[1.treated#1.post] 
 local ub_post = _b[1.treated#1.post] + invttail(e(df_r),0.025)*_se[1.treated#1.post]
 
 mat `var'_general_ci_season`s' = ((`lb_post')\(`ub_post'))  
	 
	 
*----------------------------------------------------	

*##############---------------------------


xtreg `var' ///
      treated##post1  treated##post2 i.year_season ///
			 if  year<=2015 & season==`s' & elevation>1000 & (crop_mask_v2>50|rangeland_mask_v2>50|crop_range_mask>50), fe vce(cluster mws_pixel_id_5km)
			*------------------------------------------------
   
   outreg2 using "$TABLES/DID_results_`var'_2phases_season`s'.tex", ///
   keep(1.treated#1.post1 1.treated#1.post2 total_precip_chirpsv20_s av_temp_era5_s total_precip_prev_s elev_temp temp_sq precip_sq elev_preci) replace label
  *---------------------------------------------------- 
  *---------------------------------------------------- 
  xtreg `var' ///
      treated##post1  treated##post2 i.year_season  total_precip_chirpsv20_s  av_temp_era5_s  total_precip_prev_s ///
							elev_preci elev_temp temp_sq precip_sq precip_wate    ///
			 if  year<=2015 & season==`s' & elevation>1000 & (crop_mask_v2>50|rangeland_mask_v2>50|crop_range_mask>50), fe vce(cluster mws_pixel_id_5km)
			*------------------------------------------------
   
   outreg2 using "$TABLES/DID_results_`var'_2phases_season`s'.tex", ///
   keep(1.treated#1.post1 1.treated#1.post2 total_precip_chirpsv20_s av_temp_era5_s total_precip_prev_s elev_temp temp_sq precip_sq elev_preci) append label
  *---------------------------------------------------- 
    xtreg `var' ///
      treated##post1  treated##post2 i.year_season  total_precip_chirpsv20_s  av_temp_era5_s  total_precip_prev_s ///
							elev_preci elev_temp temp_sq precip_sq precip_wate    ///
							$soil ///
			 if  year<=2015 & season==`s' & elevation>1000 & (crop_mask_v2>50|rangeland_mask_v2>50|crop_range_mask>50), fe vce(cluster mws_pixel_id_5km)
			*------------------------------------------------
   
   outreg2 using "$TABLES/DID_results_`var'_2phases_season`s'.tex", ///
   keep(1.treated#1.post1 1.treated#1.post2 total_precip_chirpsv20_s av_temp_era5_s total_precip_prev_s elev_temp temp_sq precip_sq elev_preci) append label
  *---------------------------------------------------- 
  xtreg `var' ///
      treated##post1  treated##post2 i.year_season##i.aez_09  total_precip_chirpsv20_s  av_temp_era5_s  total_precip_prev_s ///
							elev_preci elev_temp temp_sq precip_sq precip_wate    ///
							$soil ///
			 if  year<=2015 & season==`s' & elevation>1000 & (crop_mask_v2>50|rangeland_mask_v2>50|crop_range_mask>50), fe vce(cluster mws_pixel_id_5km)
			*------------------------------------------------
  
         outreg2 using "$TABLES/DID_results_`var'_2phases_season`s'.tex", ///
   keep(1.treated#1.post1 1.treated#1.post2 total_precip_chirpsv20_s av_temp_era5_s total_precip_prev_s elev_temp temp_sq precip_sq elev_preci) append label 
   
   
 * Coefficients for percentage change * 

 global scale2_`var'_1_`s' = (1/(`m'))*100 //% change formula
 global scale2_`var'_2_`s' = (1/(`m'))*100 //% change formula

  mat `var'_general_season`s' = (nullmat(`var'_general_season`s')), (_b[1.treated#1.post1]), (_b[1.treated#1.post2])

  forvalues i=1/2{
 local lb_post`i' = _b[1.treated#1.post`i'] - invttail(e(df_r),0.025)*_se[1.treated#1.post`i'] 
 local ub_post`i' = _b[1.treated#1.post`i'] + invttail(e(df_r),0.025)*_se[1.treated#1.post`i']
 }
 
 mat `var'_general_ci_season`s' =(nullmat(`var'_general_ci_season`s')), ((`lb_post1')\(`ub_post1')) , ((`lb_post2')\(`ub_post2')) 
 
 mat colnames `var'_general_season`s' = ATE ATE_post1 ATE_post2
 mat colnames `var'_general_ci_season`s' = ATE ATE_post1 ATE_post2
 

  *---------------------------------------------------- 
    test _b[1.treated#1.post1]=_b[1.treated#1.post2]
  test _b[1.treated#1.post1]=0 ,notest 
  test _b[1.treated#1.post2]=0 , accum
  
*#################################################################################################################

 *=========================================================================
 *  DROUGHT ANALYSIS 
 *=========================================================================
*============================================================
*========================================================
	   	  
	   foreach x in d_SPI_15 drought_SPI_2  {

 	xtreg `var' ///
      treated##post  treated##`x'##post treated##`x' i.year_season ///
    if  year<=2015 & season==`s' & elevation>1000 & (crop_mask_v2>50|rangeland_mask_v2>50|crop_range_mask>50), fe vce(cluster mws_pixel_id_5km)
	
	 outreg2 using "$TABLES/DID_results_`x'_`var'_season`s'.tex", keep(1.treated#1.post  1.post 1.treated#1.`x'#1.post 1.`x'#1.post 1.treated#1.`x' 1.`x' av_temp_era5_s  elev_temp temp_sq ) replace label
	
*----------------------------------------------------------------------
 	xtreg `var' ///
      treated##post  treated##`x'##post treated##`x' i.year_season  ///
	         av_temp_era5_s elev_temp temp_sq   ///   
    if  year<=2015 & season==`s' & elevation>1000 & (crop_mask_v2>50|rangeland_mask_v2>50|crop_range_mask>50), fe vce(cluster mws_pixel_id_5km)
	
	 outreg2 using "$TABLES/DID_results_`x'_`var'_season`s'.tex", keep(1.treated#1.post  1.post 1.treated#1.`x'#1.post 1.`x'#1.post 1.treated#1.`x' 1.`x' av_temp_era5_s  elev_temp temp_sq ) append label
*----------------------------------------------------------------------
 	xtreg `var' ///
      treated##post  treated##`x'##post treated##`x' i.year_season  ///
	         av_temp_era5_s elev_temp temp_sq   ///
			 $soil ///
    if  year<=2015 & season==`s' & elevation>1000 & (crop_mask_v2>50|rangeland_mask_v2>50|crop_range_mask>50), fe vce(cluster mws_pixel_id_5km)
	
	 outreg2 using "$TABLES/DID_results_`x'_`var'_season`s'.tex", keep(1.treated#1.post  1.post 1.treated#1.`x'#1.post 1.`x'#1.post 1.treated#1.`x' 1.`x' av_temp_era5_s  elev_temp temp_sq ) append label
*----------------------------------------------------------------------
	 
	  	xtreg `var' ///
      treated##post  treated##`x'##post treated##`x' i.year_season##i.aez_09  ///
	         av_temp_era5_s elev_temp temp_sq   ///   
			 $soil ///
    if  year<=2015 & season==`s' & elevation>1000 & (crop_mask_v2>50|rangeland_mask_v2>50|crop_range_mask>50), fe vce(cluster mws_pixel_id_5km)

   outreg2 using "$TABLES/DID_results_`x'_`var'_season`s'.tex", keep(1.treated#1.post  1.post 1.treated#1.`x'#1.post 1.`x'#1.post 1.treated#1.`x' 1.`x' av_temp_era5_s  elev_temp temp_sq ) append label
	*----------------------------------------------------------------------
	
 * Coefficients for percentage change *
   
      sum `var' if treated==1 & post ==0 & `x'==0 & e(sample)==1
	  local m1 = `r(mean)'
	  
	global s_`var'_`x'_post = (1/(`m1'))*100 //% change formula
    global s_`var'_`x' = (1/(`m1'))*100 //% change formula
 

   lincom 1.treated#1.`x'#1.post+1.treated#1.post
   mat `var'_`x' = (`r(estimate)')
   mat `var'_`x'_ci = (`r(lb)' \ `r(ub)')


 *----------------------------------------------------------------------
 
  mat `var'_no_`x' =  (_b[1.treated#1.post])

 local lb_post = _b[1.treated#1.post] - invttail(e(df_r),0.025)*_se[1.treated#1.post] 
 local ub_post = _b[1.treated#1.post] + invttail(e(df_r),0.025)*_se[1.treated#1.post]
 
 
 mat `var'_no_`x'_ci =((`lb_post')\(`ub_post')) 
 

*==========================================================================
*========================================================
 * Dynamic effects on drought *
*=========================================================
	xtreg `var' ///
       treated##post1 treated##`x'##post1 treated##`x'##post2 treated##`x' i.year_season ///
    if  year<=2015 & season==`s' & elevation>1000 & (crop_mask_v2>50|rangeland_mask_v2>50|crop_range_mask>50), fe vce(cluster mws_pixel_id_5km)
	
   outreg2 using "$TABLES/DID_results_`x'_`var'_2phases_season`s'.tex", ///
   keep(1.treated#1.post1 1.treated#1.post2 1.treated#1.`x'#1.post1 1.treated#1.`x'#1.post2 1.post1 1.post2  1.`x'#1.post1 1.`x'#1.post2 1.treated#1.`x' 1.`x' av_temp_era5_s  elev_temp temp_sq ) append label
   *----------------------------------------------------------------------
   	xtreg `var' ///
       treated##post1 treated##`x'##post1 treated##`x'##post2 treated##`x' i.year_season   ///
	         av_temp_era5_s elev_temp temp_sq    ///  
    if  year<=2015 & season==`s' & elevation>1000 & (crop_mask_v2>50|rangeland_mask_v2>50|crop_range_mask>50), fe vce(cluster mws_pixel_id_5km)
	
   outreg2 using "$TABLES/DID_results_`x'_`var'_2phases_season`s'.tex", ///
   keep(1.treated#1.post1 1.treated#1.post2 1.treated#1.`x'#1.post1 1.treated#1.`x'#1.post2 1.post1 1.post2  1.`x'#1.post1 1.`x'#1.post2 1.treated#1.`x' 1.`x' av_temp_era5_s  elev_temp temp_sq ) append label
   *----------------------------------------------------------------------
      	xtreg `var' ///
       treated##post1 treated##`x'##post1 treated##`x'##post2 treated##`x' i.year_season   ///
	         av_temp_era5_s elev_temp temp_sq    ///  
			 $soil ///
    if  year<=2015 & season==`s' & elevation>1000 & (crop_mask_v2>50|rangeland_mask_v2>50|crop_range_mask>50), fe vce(cluster mws_pixel_id_5km)

   outreg2 using "$TABLES/DID_results_`x'_`var'_2phases_season`s'.tex", ///
   keep(1.treated#1.post1 1.treated#1.post2 1.treated#1.`x'#1.post1 1.treated#1.`x'#1.post2 1.post1 1.post2  1.`x'#1.post1 1.`x'#1.post2 1.treated#1.`x' 1.`x' av_temp_era5_s  elev_temp temp_sq ) append label
   *----------------------------------------------------------------------
   	xtreg `var' ///
       treated##post1  treated##`x'##post1 treated##`x'##post2 treated##`x' i.year_season##i.aez_09   ///
	         av_temp_era5_s elev_temp temp_sq    ///  
			 $soil ///
    if  year<=2015 & season==`s' & elevation>1000 & (crop_mask_v2>50|rangeland_mask_v2>50|crop_range_mask>50), fe vce(cluster mws_pixel_id_5km)
	    
		outreg2 using "$TABLES/DID_results_`x'_`var'_2phases_season`s'.tex", ///
   keep(1.treated#1.post1 1.treated#1.post2 1.treated#1.`x'#1.post1 1.treated#1.`x'#1.post2 1.post1 1.post2  1.`x'#1.post1 1.`x'#1.post2 1.treated#1.`x' 1.`x' av_temp_era5_s  elev_temp temp_sq ) append label
   *----------------------------------------------------------------------

	  sum `var' if treated==1 & post==0  & e(sample)==1
	  local m1 = `r(mean)'
	  
 global s_`var'_pos1_d = (1/(`m1'))*100 //% change formula
 global s_`var'_pos2_d = (1/(`m1'))*100 //% change formula
 global s_`var'_pos2_`x' = (1/(`m1'))*100 //% change formula
 global s_`var'_pos1_`x' = (1/(`m1'))*100 //% change formula
 
/*
  lincom 1.treated#1.`x'#1.post1+1.treated#1.post1 //cannot identify severe drought in the post1 period
   mat `var'_`x' = (nullmat(`var'_`x')), (`r(estimate)')
   mat `var'_`x'_ci = (nullmat(`var'_`x'_ci)), (`r(lb)' \ `r(ub)')
*/
  
 lincom 1.treated#1.`x'#1.post2+1.treated#1.post2
   mat `var'_`x' = (nullmat(`var'_`x')), (.), (`r(estimate)')
   mat `var'_`x'_ci = (nullmat(`var'_`x'_ci)), ((.)\(.)), (`r(lb)' \ `r(ub)')
mat colnames `var'_`x' = ATE ATE_post1 ATE_post2
mat colnames `var'_`x'_ci = ATE ATE_post1 ATE_post2


  mat `var'_no_`x' = (nullmat(`var'_no_`x')), (_b[1.treated#1.post1]), (_b[1.treated#1.post2])

  forvalues i=1/2{
 local lb_post`i' = _b[1.treated#1.post`i'] - invttail(e(df_r),0.025)*_se[1.treated#1.post`i'] 
 local ub_post`i' = _b[1.treated#1.post`i'] + invttail(e(df_r),0.025)*_se[1.treated#1.post`i']
 }
 
 mat `var'_no_`x'_ci =(nullmat(`var'_no_`x'_ci)), ((`lb_post1')\(`ub_post1')) , ((`lb_post2')\(`ub_post2')) 
 
 mat colnames `var'_no_`x' = ATE ATE_post1 ATE_post2
 mat colnames `var'_no_`x'_ci = ATE ATE_post1 ATE_post2
 
 *------------------------------------------------------------------------
 *------------------------------------------------------------------------
   test _b[1.treated#1.post1]=_b[1.treated#1.post2]
  test _b[1.treated#1.post1]=0 ,notest 
  test _b[1.treated#1.post2]=0 , accum
 *-----------------------------------------------------
 
 }
}
 }
 
    	 forvalues s =3/3{
 foreach var of varlist  evi_5km_s sif_s gosif_gpp_s {
    
	   foreach x in     drought_SPI_1 {

 	xtreg `var' ///
      treated##post  treated##`x'##post treated##`x' i.year_season ///
    if  year<=2015 & season==`s' & elevation>1000 & (crop_mask_v2>50|rangeland_mask_v2>50|crop_range_mask>50), fe vce(cluster mws_pixel_id_5km)
	
	 outreg2 using "$TABLES/DID_results_`x'_`var'_season`s'.tex", keep(1.treated#1.post  1.post 1.treated#1.`x'#1.post 1.`x'#1.post 1.treated#1.`x' 1.`x' av_temp_era5_s  elev_temp temp_sq ) replace label
	
*----------------------------------------------------------------------
 	xtreg `var' ///
      treated##post  treated##`x'##post treated##`x' i.year_season  ///
	         av_temp_era5_s elev_temp temp_sq   ///   
    if  year<=2015 & season==`s' & elevation>1000 & (crop_mask_v2>50|rangeland_mask_v2>50|crop_range_mask>50), fe vce(cluster mws_pixel_id_5km)
	
	 outreg2 using "$TABLES/DID_results_`x'_`var'_season`s'.tex", keep(1.treated#1.post  1.post 1.treated#1.`x'#1.post 1.`x'#1.post 1.treated#1.`x' 1.`x' av_temp_era5_s  elev_temp temp_sq ) append label
*----------------------------------------------------------------------
 	xtreg `var' ///
      treated##post  treated##`x'##post treated##`x' i.year_season  ///
	         av_temp_era5_s elev_temp temp_sq   ///
			 $soil ///
    if  year<=2015 & season==`s' & elevation>1000 & (crop_mask_v2>50|rangeland_mask_v2>50|crop_range_mask>50), fe vce(cluster mws_pixel_id_5km)
	
	 outreg2 using "$TABLES/DID_results_`x'_`var'_season`s'.tex", keep(1.treated#1.post  1.post 1.treated#1.`x'#1.post 1.`x'#1.post 1.treated#1.`x' 1.`x' av_temp_era5_s  elev_temp temp_sq ) append label
*----------------------------------------------------------------------
	 
	  	xtreg `var' ///
      treated##post  treated##`x'##post treated##`x' i.year_season##i.aez_09  ///
	         av_temp_era5_s elev_temp temp_sq   ///   
			 $soil ///
    if  year<=2015 & season==`s' & elevation>1000 & (crop_mask_v2>50|rangeland_mask_v2>50|crop_range_mask>50), fe vce(cluster mws_pixel_id_5km)

   outreg2 using "$TABLES/DID_results_`x'_`var'_season`s'.tex", keep(1.treated#1.post  1.post 1.treated#1.`x'#1.post 1.`x'#1.post 1.treated#1.`x' 1.`x' av_temp_era5_s  elev_temp temp_sq ) append label
	*----------------------------------------------------------------------
	
 * Coefficients for percentage change *
   
      sum `var' if treated==1 & post ==0 & `x'==0 & e(sample)==1
	  local m1 = `r(mean)'
	  
	global s_`var'_`x'_post = (1/(`m1'))*100 //% change formula
    global s_`var'_`x' = (1/(`m1'))*100 //% change formula
 

   lincom 1.treated#1.`x'#1.post+1.treated#1.post
   mat `var'_`x' = (`r(estimate)')
   mat `var'_`x'_ci = (`r(lb)' \ `r(ub)')


 *----------------------------------------------------------------------
 
  mat `var'_no_`x' =  (_b[1.treated#1.post])

 local lb_post = _b[1.treated#1.post] - invttail(e(df_r),0.025)*_se[1.treated#1.post] 
 local ub_post = _b[1.treated#1.post] + invttail(e(df_r),0.025)*_se[1.treated#1.post]
 
 
 mat `var'_no_`x'_ci =((`lb_post')\(`ub_post')) 
 

*==========================================================================
*========================================================
 * Dynamic effects on drought *
*=========================================================
	xtreg `var' ///
       treated##post1 treated##`x'##post1 treated##`x'##post2 treated##`x' i.year_season ///
    if  year<=2015 & season==`s' & elevation>1000 & (crop_mask_v2>50|rangeland_mask_v2>50|crop_range_mask>50), fe vce(cluster mws_pixel_id_5km)
	
   outreg2 using "$TABLES/DID_results_`x'_`var'_2phases_season`s'.tex", ///
   keep(1.treated#1.post1 1.treated#1.post2 1.treated#1.`x'#1.post1 1.treated#1.`x'#1.post2 1.post1 1.post2  1.`x'#1.post1 1.`x'#1.post2 1.treated#1.`x' 1.`x' av_temp_era5_s  elev_temp temp_sq ) append label
   *----------------------------------------------------------------------
   	xtreg `var' ///
       treated##post1 treated##`x'##post1 treated##`x'##post2 treated##`x' i.year_season   ///
	         av_temp_era5_s elev_temp temp_sq    ///  
    if  year<=2015 & season==`s' & elevation>1000 & (crop_mask_v2>50|rangeland_mask_v2>50|crop_range_mask>50), fe vce(cluster mws_pixel_id_5km)
	
   outreg2 using "$TABLES/DID_results_`x'_`var'_2phases_season`s'.tex", ///
   keep(1.treated#1.post1 1.treated#1.post2 1.treated#1.`x'#1.post1 1.treated#1.`x'#1.post2 1.post1 1.post2  1.`x'#1.post1 1.`x'#1.post2 1.treated#1.`x' 1.`x' av_temp_era5_s  elev_temp temp_sq ) append label
   *----------------------------------------------------------------------
      	xtreg `var' ///
       treated##post1 treated##`x'##post1 treated##`x'##post2 treated##`x' i.year_season   ///
	         av_temp_era5_s elev_temp temp_sq    ///  
			 $soil ///
    if  year<=2015 & season==`s' & elevation>1000 & (crop_mask_v2>50|rangeland_mask_v2>50|crop_range_mask>50), fe vce(cluster mws_pixel_id_5km)

   outreg2 using "$TABLES/DID_results_`x'_`var'_2phases_season`s'.tex", ///
   keep(1.treated#1.post1 1.treated#1.post2 1.treated#1.`x'#1.post1 1.treated#1.`x'#1.post2 1.post1 1.post2  1.`x'#1.post1 1.`x'#1.post2 1.treated#1.`x' 1.`x' av_temp_era5_s  elev_temp temp_sq ) append label
   *----------------------------------------------------------------------
   	xtreg `var' ///
       treated##post1  treated##`x'##post1 treated##`x'##post2 treated##`x' i.year_season##i.aez_09   ///
	         av_temp_era5_s elev_temp temp_sq    ///  
			 $soil ///
    if  year<=2015 & season==`s' & elevation>1000 & (crop_mask_v2>50|rangeland_mask_v2>50|crop_range_mask>50), fe vce(cluster mws_pixel_id_5km)
	    
		outreg2 using "$TABLES/DID_results_`x'_`var'_2phases_season`s'.tex", ///
   keep(1.treated#1.post1 1.treated#1.post2 1.treated#1.`x'#1.post1 1.treated#1.`x'#1.post2 1.post1 1.post2  1.`x'#1.post1 1.`x'#1.post2 1.treated#1.`x' 1.`x' av_temp_era5_s  elev_temp temp_sq ) append label
   *----------------------------------------------------------------------

	  sum `var' if treated==1 & post==0  & e(sample)==1
	  local m1 = `r(mean)'
	  
 global s_`var'_pos1_d = (1/(`m1'))*100 //% change formula
 global s_`var'_pos2_d = (1/(`m1'))*100 //% change formula
 global s_`var'_pos2_`x' = (1/(`m1'))*100 //% change formula
 global s_`var'_pos1_`x' = (1/(`m1'))*100 //% change formula
 

  lincom 1.treated#1.`x'#1.post1+1.treated#1.post1
   mat `var'_`x' = (nullmat(`var'_`x')), (`r(estimate)')
   mat `var'_`x'_ci = (nullmat(`var'_`x'_ci)), (`r(lb)' \ `r(ub)')

  
 lincom 1.treated#1.`x'#1.post2+1.treated#1.post2
   mat `var'_`x' = (nullmat(`var'_`x')), (`r(estimate)')
   mat `var'_`x'_ci = (nullmat(`var'_`x'_ci)), (`r(lb)' \ `r(ub)')
mat colnames `var'_`x' = ATE ATE_post1 ATE_post2
mat colnames `var'_`x'_ci = ATE ATE_post1 ATE_post2


  mat `var'_no_`x' = (nullmat(`var'_no_`x')), (_b[1.treated#1.post1]), (_b[1.treated#1.post2])

  forvalues i=1/2{
 local lb_post`i' = _b[1.treated#1.post`i'] - invttail(e(df_r),0.025)*_se[1.treated#1.post`i'] 
 local ub_post`i' = _b[1.treated#1.post`i'] + invttail(e(df_r),0.025)*_se[1.treated#1.post`i']
 }
 
 mat `var'_no_`x'_ci =(nullmat(`var'_no_`x'_ci)), ((`lb_post1')\(`ub_post1')) , ((`lb_post2')\(`ub_post2')) 
 
 mat colnames `var'_no_`x' = ATE ATE_post1 ATE_post2
 mat colnames `var'_no_`x'_ci = ATE ATE_post1 ATE_post2
 
 *------------------------------------------------------------------------
 *------------------------------------------------------------------------
   test _b[1.treated#1.post1]=_b[1.treated#1.post2]
  test _b[1.treated#1.post1]=0 ,notest 
  test _b[1.treated#1.post2]=0 , accum
 *-----------------------------------------------------
 
 }
}
 }
 
****************************
* Robustness check - alternative landcover map 

 forvalues s =3/3{
 foreach var of varlist  evi_5km_s sif_s gosif_gpp_s {
   
	   xtset mws_pixel_id_5km year_season 
  
  *###--------------------------------
	
xtreg `var' ///
      treated##post  i.year_season ///
			 if  year<=2015 & season==`s' & elevation>1000 & (crop_frac>0.50|range_frac>0.50|crop_range_frac>0.50), fe vce(cluster mws_pixel_id_5km)
			*------------------------------------------------
	
	outreg2 using "$TABLES/DID_results_`var'_lcover2_season`s'.tex", ///
	 keep(1.treated#1.post total_precip_chirpsv20_s av_temp_era5_s total_precip_prev_s elev_temp temp_sq precip_sq elev_preci) replace label
*----------------------------------------------------	 
xtreg `var' ///
      treated##post  i.year_season   total_precip_chirpsv20_s  av_temp_era5_s total_precip_prev_s ///
							elev_preci elev_temp temp_sq precip_sq precip_wate   ///
			 if  year<=2015 & season==`s' & elevation>1000 & (crop_frac>0.50|range_frac>0.50|crop_range_frac>0.50), fe vce(cluster mws_pixel_id_5km)
			*------------------------------------------------
	
	outreg2 using "$TABLES/DID_results_`var'_lcover2_season`s'.tex", ///
	 keep(1.treated#1.post total_precip_chirpsv20_s av_temp_era5_s total_precip_prev_s elev_temp temp_sq precip_sq elev_preci) append label
*----------------------------------------------------	
xtreg `var' ///
      treated##post  i.year_season   total_precip_chirpsv20_s  av_temp_era5_s total_precip_prev_s ///
							elev_preci elev_temp temp_sq precip_sq precip_wate   ///
							$soil ///
			 if  year<=2015 & season==`s' & elevation>1000 & (crop_frac>0.50|range_frac>0.50|crop_range_frac>0.50), fe vce(cluster mws_pixel_id_5km)
			*------------------------------------------------
	
	outreg2 using "$TABLES/DID_results_`var'_lcover2_season`s'.tex", ///
	 keep(1.treated#1.post total_precip_chirpsv20_s av_temp_era5_s total_precip_prev_s elev_temp temp_sq precip_sq elev_preci) append label
*----------------------------------------------------	
xtreg `var' ///
      treated##post  i.year_season##i.aez_09   total_precip_chirpsv20_s  av_temp_era5_s total_precip_prev_s ///
							elev_preci elev_temp temp_sq precip_sq precip_wate   ///
							$soil ///
			 if  year<=2015 & season==`s' & elevation>1000 & (crop_frac>0.50|range_frac>0.50|crop_range_frac>0.50), fe vce(cluster mws_pixel_id_5km)
			*------------------------------------------------
  	outreg2 using "$TABLES/DID_results_`var'_lcover2_season`s'.tex", ///
	 keep(1.treated#1.post total_precip_chirpsv20_s av_temp_era5_s total_precip_prev_s elev_temp temp_sq precip_sq elev_preci) append label ///
	 addnote(treatment group mean in the post period = `m')

*##############---------------------------


xtreg `var' ///
      treated##post1  treated##post2 i.year_season ///
			 if  year<=2015 & season==`s' & elevation>1000 & (crop_frac>0.50|range_frac>0.50|crop_range_frac>0.50), fe vce(cluster mws_pixel_id_5km)
			*------------------------------------------------
   
   outreg2 using "$TABLES/DID_results_`var'_2phases_lcover2_season`s'.tex", ///
   keep(1.treated#1.post1 1.treated#1.post2 total_precip_chirpsv20_s av_temp_era5_s total_precip_prev_s elev_temp temp_sq precip_sq elev_preci) replace label
  *---------------------------------------------------- 
  *---------------------------------------------------- 
  xtreg `var' ///
      treated##post1  treated##post2 i.year_season  total_precip_chirpsv20_s  av_temp_era5_s  total_precip_prev_s ///
							elev_preci elev_temp temp_sq precip_sq precip_wate    ///
			 if  year<=2015 & season==`s' & elevation>1000 & (crop_frac>0.50|range_frac>0.50|crop_range_frac>0.50), fe vce(cluster mws_pixel_id_5km)
			*------------------------------------------------
   
   outreg2 using "$TABLES/DID_results_`var'_2phases_lcover2_season`s'.tex", ///
   keep(1.treated#1.post1 1.treated#1.post2 total_precip_chirpsv20_s av_temp_era5_s total_precip_prev_s elev_temp temp_sq precip_sq elev_preci) append label
  *---------------------------------------------------- 
    xtreg `var' ///
      treated##post1  treated##post2 i.year_season  total_precip_chirpsv20_s  av_temp_era5_s  total_precip_prev_s ///
							elev_preci elev_temp temp_sq precip_sq precip_wate    ///
							$soil ///
			 if  year<=2015 & season==`s' & elevation>1000 & (crop_frac>0.50|range_frac>0.50|crop_range_frac>0.50), fe vce(cluster mws_pixel_id_5km)
			*------------------------------------------------
   
   outreg2 using "$TABLES/DID_results_`var'_2phases_lcover2_season`s'.tex", ///
   keep(1.treated#1.post1 1.treated#1.post2 total_precip_chirpsv20_s av_temp_era5_s total_precip_prev_s elev_temp temp_sq precip_sq elev_preci) append label
  *---------------------------------------------------- 
  xtreg `var' ///
      treated##post1  treated##post2 i.year_season##i.aez_09  total_precip_chirpsv20_s  av_temp_era5_s  total_precip_prev_s ///
							elev_preci elev_temp temp_sq precip_sq precip_wate    ///
							$soil ///
			 if  year<=2015 & season==`s' & elevation>1000 & (crop_frac>0.50|range_frac>0.50|crop_range_frac>0.50), fe vce(cluster mws_pixel_id_5km)
			*------------------------------------------------
  
         outreg2 using "$TABLES/DID_results_`var'_2phases_lcover2_season`s'.tex", ///
   keep(1.treated#1.post1 1.treated#1.post2 total_precip_chirpsv20_s av_temp_era5_s total_precip_prev_s elev_temp temp_sq precip_sq elev_preci) append label 
   
  
*#################################################################################################################

 *=========================================================================
 *  DROUGHT ANALYSIS 
 *=========================================================================
*============================================================
*========================================================
	   	  
	   foreach x in drought_SPI_2  {

 	xtreg `var' ///
      treated##post  treated##`x'##post treated##`x' i.year_season ///
    if  year<=2015 & season==`s' & elevation>1000 & (crop_frac>0.50|range_frac>0.50|crop_range_frac>0.50), fe vce(cluster mws_pixel_id_5km)
	
	 outreg2 using "$TABLES/DID_results_`x'_`var'_lcover2_season`s'.tex", keep(1.treated#1.post  1.post 1.treated#1.`x'#1.post 1.`x'#1.post 1.treated#1.`x' 1.`x' av_temp_era5_s  elev_temp temp_sq ) replace label
	
*----------------------------------------------------------------------
 	xtreg `var' ///
      treated##post  treated##`x'##post treated##`x' i.year_season  ///
	         av_temp_era5_s elev_temp temp_sq   ///   
    if  year<=2015 & season==`s' & elevation>1000 & (crop_frac>0.50|range_frac>0.50|crop_range_frac>0.50), fe vce(cluster mws_pixel_id_5km)
	
	 outreg2 using "$TABLES/DID_results_`x'_`var'_lcover2_season`s'.tex", keep(1.treated#1.post  1.post 1.treated#1.`x'#1.post 1.`x'#1.post 1.treated#1.`x' 1.`x' av_temp_era5_s  elev_temp temp_sq ) append label
*----------------------------------------------------------------------
 	xtreg `var' ///
      treated##post  treated##`x'##post treated##`x' i.year_season  ///
	         av_temp_era5_s elev_temp temp_sq   ///
			 $soil ///
    if  year<=2015 & season==`s' & elevation>1000 & (crop_frac>0.50|range_frac>0.50|crop_range_frac>0.50), fe vce(cluster mws_pixel_id_5km)
	
	 outreg2 using "$TABLES/DID_results_`x'_`var'_lcover2_season`s'.tex", keep(1.treated#1.post  1.post 1.treated#1.`x'#1.post 1.`x'#1.post 1.treated#1.`x' 1.`x' av_temp_era5_s  elev_temp temp_sq ) append label
*----------------------------------------------------------------------
	 
	  	xtreg `var' ///
      treated##post  treated##`x'##post treated##`x' i.year_season##i.aez_09  ///
	         av_temp_era5_s elev_temp temp_sq   ///   
			 $soil ///
    if  year<=2015 & season==`s' & elevation>1000 & (crop_frac>0.50|range_frac>0.50|crop_range_frac>0.50), fe vce(cluster mws_pixel_id_5km)

   outreg2 using "$TABLES/DID_results_`x'_`var'_lcover2_season`s'.tex", keep(1.treated#1.post  1.post 1.treated#1.`x'#1.post 1.`x'#1.post 1.treated#1.`x' 1.`x' av_temp_era5_s  elev_temp temp_sq ) append label
	*----------------------------------------------------------------------

*==========================================================================
*========================================================
 * Dynamic effects on drought *
*=========================================================
	xtreg `var' ///
       treated##post1 treated##`x'##post1 treated##`x'##post2 treated##`x' i.year_season ///
    if  year<=2015 & season==`s' & elevation>1000 & (crop_frac>0.50|range_frac>0.50|crop_range_frac>0.50), fe vce(cluster mws_pixel_id_5km)
	
   outreg2 using "$TABLES/DID_results_`x'_`var'_2phases_lcover2_season`s'.tex", ///
   keep(1.treated#1.post1 1.treated#1.post2 1.treated#1.`x'#1.post1 1.treated#1.`x'#1.post2 1.post1 1.post2  1.`x'#1.post1 1.`x'#1.post2 1.treated#1.`x' 1.`x' av_temp_era5_s  elev_temp temp_sq ) append label
   *----------------------------------------------------------------------
   	xtreg `var' ///
       treated##post1 treated##`x'##post1 treated##`x'##post2 treated##`x' i.year_season   ///
	         av_temp_era5_s elev_temp temp_sq    ///  
    if  year<=2015 & season==`s' & elevation>1000 & (crop_frac>0.50|range_frac>0.50|crop_range_frac>0.50), fe vce(cluster mws_pixel_id_5km)
	
   outreg2 using "$TABLES/DID_results_`x'_`var'_2phases_lcover2_season`s'.tex", ///
   keep(1.treated#1.post1 1.treated#1.post2 1.treated#1.`x'#1.post1 1.treated#1.`x'#1.post2 1.post1 1.post2  1.`x'#1.post1 1.`x'#1.post2 1.treated#1.`x' 1.`x' av_temp_era5_s  elev_temp temp_sq ) append label
   *----------------------------------------------------------------------
      	xtreg `var' ///
       treated##post1 treated##`x'##post1 treated##`x'##post2 treated##`x' i.year_season   ///
	         av_temp_era5_s elev_temp temp_sq    ///  
			 $soil ///
    if  year<=2015 & season==`s' & elevation>1000 & (crop_frac>0.50|range_frac>0.50|crop_range_frac>0.50), fe vce(cluster mws_pixel_id_5km)

   outreg2 using "$TABLES/DID_results_`x'_`var'_2phases_lcover2_season`s'.tex", ///
   keep(1.treated#1.post1 1.treated#1.post2 1.treated#1.`x'#1.post1 1.treated#1.`x'#1.post2 1.post1 1.post2  1.`x'#1.post1 1.`x'#1.post2 1.treated#1.`x' 1.`x' av_temp_era5_s  elev_temp temp_sq ) append label
   *----------------------------------------------------------------------
   	xtreg `var' ///
       treated##post1  treated##`x'##post1 treated##`x'##post2 treated##`x' i.year_season##i.aez_09   ///
	         av_temp_era5_s elev_temp temp_sq    ///  
			 $soil ///
    if  year<=2015 & season==`s' & elevation>1000 & (crop_frac>0.50|range_frac>0.50|crop_range_frac>0.50), fe vce(cluster mws_pixel_id_5km)
	    
		outreg2 using "$TABLES/DID_results_`x'_`var'_2phases_lcover2_season`s'.tex", ///
   keep(1.treated#1.post1 1.treated#1.post2 1.treated#1.`x'#1.post1 1.treated#1.`x'#1.post2 1.post1 1.post2  1.`x'#1.post1 1.`x'#1.post2 1.treated#1.`x' 1.`x' av_temp_era5_s  elev_temp temp_sq ) append label
   *----------------------------------------------------------------------
 
 }
}
 }
 
    	 
*************************************************************
* PLOTTING RESULTS GRAPHICALLY *
*************************************************************	 
*----------------------------------------------
 
	coefplot (matrix(evi_general_season3), ci(evi_general_ci_season3) rescale(ATE=$scale1_EVI_3 ATE_post1=$scale2_EVI_1_3 ATE_post2=$scale2_EVI_2_3)  ///
            mlw(medthick) mlc(orange_red) mc(gray) mcolor(orange_red)  ciopts(recast(rcap) lwidth(medthick) lcolor(orange_red))  msymbol(Dh) msize(medlarge) offset(0.10) label(EVI)) ///
         (matrix(gosif_gpp_s_general_season3), ci(gosif_gpp_s_general_ci_season3) rescale(ATE=$scale1_gosif_gpp_s_3 ATE_post1=$scale2_gosif_gpp_s_1_3 ATE_post2=$scale2_gosif_gpp_s_2_3)  ///
		  mlw(medthick) mlc(black) mcolor(black)  ciopts(recast(rcap) lwidth(medthick) lcolor(black)) msymbol(Oh) msize(medlarge)  offset(-0.10) label(GOSIF-GPP)) ///
			 , bylabel({bf:A.} Total effects) ///
 ||  (matrix(evi_no_drought), ci(evi_no_drought_ci) rescale(ATE=$s_no_drought_SPI_1_EVI ATE_post1=$scale_post1_EVI_d ATE_post2=$scale_post2_EVI_d) ) ///
      (matrix(gosif_gpp_s_no_drought_SPI_1), ci(gosif_gpp_s_no_drought_SPI_1_ci) rescale(ATE=$s_gosif_gpp_s_drought_SPI_1 ATE_post1=$s_gosif_gpp_s_pos1_d ATE_post2=$s_gosif_gpp_s_pos2_d) ) ///
		 , bylabel({bf:B.} No drought areas ) /// 
 ||  (matrix(evi_drought_SPI_2), ci(evi_drought_SPI_2_ci) rescale(ATE=$s_drought_SPI_2_post_EVI  ATE_post2=$s_drought_SPI_2_post2_EVI) ) ///
        (matrix(gosif_gpp_s_drought_SPI_2), ci(gosif_gpp_s_drought_SPI_2_ci) rescale(ATE=$s_gosif_gpp_s_drought_SPI_2_post  ATE_post2=$s_gosif_gpp_s_pos2_drought_SPI_2) ) ///  
			 , bylabel({bf:C.} SPI < -2 ) ///
    ||  (matrix(evi_drought_SPI_15), ci(evi_drought_SPI_15_ci) rescale(ATE=$s_drought_SPI_15_post_EVI  ATE_post2=$s_drought_SPI_15_post2_EVI) ) ///
        (matrix(gosif_gpp_s_d_SPI_15), ci(gosif_gpp_s_d_SPI_15_ci) rescale(ATE=$s_gosif_gpp_s_d_SPI_15_post  ATE_post2=$s_gosif_gpp_s_pos2_d_SPI_15) ) ///  
			 , bylabel({bf:D.} SPI < -1.5 ) ///
   || (matrix(evi_drought_SPI_1), ci(evi_drought_SPI_1_ci) rescale(ATE=$s_drought_SPI_1_post_EVI  ATE_post1=$s_drought_SPI_1_post1_EVI ATE_post2=$s_drought_SPI_1_post2_EVI) )  ///
      (matrix(gosif_gpp_s_drought_SPI_1), ci(gosif_gpp_s_drought_SPI_1_ci) rescale(ATE=$s_gosif_gpp_s_drought_SPI_1_post ATE_post1=$s_gosif_gpp_s_pos1_drought_SPI_1 ATE_post2=$s_gosif_gpp_s_pos2_drought_SPI_1) )  ///
		 , bylabel({bf:E.} SPI < -1 ) ///
	|| 		, byopts(im(vsmall) xrescale cols(5))  mlabel mlabp(2) mlabgap(2) format(%9.2f) mlabsize(medium) ///
	 xline(0, lcolor(gray%70) lpattern(shortdash) lwidth(medthick)) ///
	 xlabel(0 "0" 5 "5"  15 "15" 25 "25", labsize(medium) format(%01.0f)) legend(pos(6) rows(1)  size(medium)) ///
	 coeflabels(ATE = "% change (2011-2015)" ATE_post1="% change (2011-2013)" ATE_post2="% change (2014-2015)", labsize(small)) ///
	 ylabel(, labsize(medlarge)) 
	 
		 graph export "$GRAPHS/Fig4_Main_results_all_drought_cutoffs.tif", replace

*---------------------------------------------------------------	 
	 
*==================================================
	 

